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ABSTRACT 


A  procedure  is  presented  to  perform  a  contact  analysis  of  spiral  bevel  gears  in  order  to 
predict  the  contact  path  and  the  load  sharing  as  the  gears  roll  through  mesh.  The  approach 
utilizes  recent  advances  in  automated  contact  methods  for  nonlinear  finite  element  analysis. 
A  sector  of  the  pinion  and  gear  is  modeled  consisting  of  three  pinion  teeth  and  four  gear 
teeth  in  mesh.  Calculation  of  the  contact  force  and  stresses  through  the  gear  meshing  cycle 
are  demonstrated.  Summary  of  the  results  are  presented  using  3-dimensional  plots  and 
tables.  Issues  relating  to  solution  convergence  and  requirements  for  running  large  finite 
element  analysis  on  a  supercomputer  are  discussed. 


CHAPTER  I 


INTRODUCTION 

Spiral  bevel  gears  are  important  elements  for  transmitting  power  between  intersecting  shafts. 
They  are  commonly  used  in  applications  that  require  high  load  capacity  at  high  operating  speeds. 
One  such  application  is  in  helicopter  transmission  systems.  Aircraft  designers  are  continually 
required  to  improve  performance.  Reduced  weight,  size,  and  cost  with  increased  power, 
efficiency,  and  reliability  are  constantly  being  sought. 

Prior  research  has  focussed  on  various  aspects  of  spiral  bevel  gear  operation.  Much  has 
been  done  on  spiral  bevel  gear  geometry  to  reduce  noise  and  vibration,  kinematic  error,  improve 
manufacturability,  and  inspection.  Stress  analysis  is  another  important  area  of  ongoing 
research.  Accurate  prediction  of  contact  stresses  and  tooth  root/fillet  stresses  are  important  to 
increase  reliability  and  reduce  weight. 

Much  effort  has  focussed  on  predicting  stresses  in  gears  with  the  finite  element  method. 
Most  of  this  work  has  involved  parallel  axis  gears  with  two  dimensional  models.  Only  a  few 
researchers  have  investigated  finite  element  analysis  of  spiral  bevel  gears  (ref.  1-4).  In  reference 
4,  finite  element  analysis  was  done  on  a  single  spiral  bevel  gear  tooth  using  an  assumed  contact 
stress  distribution.  In  that  model,  contact  stresses  were  not  evaluated. 

For  parallel  axis  gears,  a  closed  form  solution  exists  which  determines  the  surface 
coordinates  of  a  tooth.  This  is  then  used  as  input  to  finite  element  methods.  For  spiral  bevel 
gears  there  is  no  closed  form  solution.  Therefore,  the  kinematics  of  the  cutting  or  grinding 
machinery  is  utilized  to  numerically  describe  the  surface  coordinates  of  the  gear  tooth. 

The  research  reported  herein  is  based  on  the  extension  of  work  presented  in  references  4-7. 
A  model  that  has  three  pinion  teeth  and  four  gear  teeth  has  been  developed  based  on  gear 
manufacturing  kinematics  for  a  single  tooth  on  each  the  pinion  and  the  gear. 
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The  objective  of  this  research  is  to  study  the  contact  path  and  load  sharing  in  these  gears 
when  contact  occurs  on  multiple  teeth  in  mesh.  This  is  done  by  performing  a  static  analysis  at 
different  incremental  rotations.  A  nonlinear  approach  is  required  due  to  large  displacements 
associated  with  gear  rotation  and  nonlinear  boundary  conditions  associated  with  the  gear  tooth 
surface  contact.  Also  evaluated  are  the  contact  stresses  and  fillet  stresses. 
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CHAPTER  n 


GEAR  SURFACE  GEOMETRY 


Briefly  described  in  this  chapter  are  the  gear  manufacturing  process,  the  kinematics  of  the 
manufacturing  process,  tooth  surface  coordinates  solution  technique,  surface  rotations  of  the  gear 
and  pinion,  and  different  orientations  required  for  the  spiral  bevel  gears  to  mesh  with  each  other. 

2.1  Gear  Manufacture 

Machinery  for  the  manufacture  of  spiral  bevel  gears  is  provided  by  the  Gleason  Works, 
Rochester,  NY.  These  machines  are  preferred  over  gear  bobbing  machines  because  they  can 
be  used  for  both  milling  and  grinding  operations.  Grinding  is  important  in  producing  hardened 
high  quality  aircraft  gears. 

This  machine  consists  mainly  of  three  parts:  the  machine  frame,  the  cradle,  and  the  sliding 
base  (ref.  8).  The  cradle  with  the  head  cutter  mounted  on  it,  slowly  rotates  about  its  axis,  as 
does  the  gear  which  is  being  cut.  During  this  motion  the  gear  surface  is  generated.  As  the 
cutter  disengages  from  the  workpiece,  the  cradle  reverses  rapidly  and  the  sliding  base  translates 
with  respect  to  the  cutter  to  index  the  workpiece  for  the  next  cutting  cycle.  This  sequence  is 
repeated  until  the  last  tooth  is  cut. 

The  head  cutter  is  mounted  on  the  cradle  with  an  offset  from  the  cradle  center.  This  allows 
adjustment  of  the  axial  distance  between  the  cutter  center  and  the  machine  center.  The 
adjustment  of  the  angular  position  between  the  two  axes  provides  the  desired  spiral  angle.  The 
shape  of  the  blades  of  the  head  cutter  are  typically  straight  lines  that  rotate  about  their  own  axis 
at  a  speed  for  efficient  metal  removal.  The  rotational  speed  of  the  cutting  head  is  independent 
of  the  cradle  or  workpiece  rotations  (ref.  9). 
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The  pinion  is  typically  cut  one  side  at  a  time,  whereas  the  gear  is  cut  both  sides 
simultaneously.  Spiral  bevel  gears  can  be  either  left  or  right  handed.  In  a  left  handed  spiral 
bevel  gear,  the  tooth  spirals  to  the  left  while  looking  from  the  apex  of  cone  towards  the  gear. 
Whereas,  for  a  right  handed  spiral  bevel  gear  the  tooth  spirals  to  the  right. 

2.2  Tooth  Surface  Coordinate  Equations 

The  system  of  equations,  required  to  define  the  coordinates  of  spiral  bevel  gear  tooth 
surfaces,  were  derived  in  reference  4,  and  are  briefly  summarized  here. 

The  first  equation  is  the  equation  of  meshing.  This  equation  is  based  on  the  kinematics  of 
manufacture  and  the  machine  tool  settings.  The  equation  of  meshing  requires  that  the  relative 
velocity  between  a  point  on  the  cutting  surface  and  the  same  point  on  the  pinion  being  cut  must 
be  perpendicular  to  the  cutting  surface  normal  (ref.  9). 

n  •  V  =  0 

where,  n  is  the  normal  vector  from  the  cutter  surface  and 

V  is  the  relative  velocity  between  the  cutter  and  the  workpiece  surfaces 
at  the  specified  location. 

This  equation  is  developed  in  terms  of  machine  tool  coordinates  U,  0  and  (t>^ 

cotij;  -  U  cosijf 
U  sinijf  sin0 
U  sinil;  COS0 
1 

where,  U  is  the  generating  cone  surface  coordinate  used  to  locate  a  point  along  the  length 

of  the  cutting  head 

0  is  the  generating  cone  surface  coordinate  used  for  rotational  orientation  of  a 
point  on  the  cutting  head 

is  the  rotated  orientation  of  the  cutter  as  it  swings  on  the  cradle 
r  is  the  radius  of  the  generating  cone  surface  and  is  described  by  the  following 
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equations  (ref.  4): 

Since  the  kinematic  motion  of  cutting  a  gear  is  equivalent  to  the  cutting  head  meshing  with  a 
simulated  crown  gear,  an  equation  of  meshing  can  be  written  in  terms  of  a  point  on  the  cutting 
head  i.e  in  terms  of  U,  0  and  4.  The  equation  of  meshing  for  straight  sided  cutters  with  a 
constant  ratio  of  roll  between  the  cutter  and  the  workpiece  is  given  by  (ref.  1): 

(U  -  r  cotijf  cosij;)  cosy  sinr 
+  -  siny)  cosilr  sin0 

T  cosy  siny  sin(g  -  4>£,)  .2) 

±  (cosy  siny  +  siny  cosy  cost)  ^  ' 

-  siny  cosijr  sint  =  0 

The  upper  and  lower  signs  are  for  left  and  right  hand  gears  respectively.  The  following 

machine  tool  settings  are  defined: 

T  (0  +  q±(/)J 

q  cradle  angle 

7  root  angle  of  workpiece 

Eta  machining  offset 

Lta  vector  sum  of  change  of  machine  center  to  back  and  the  sliding  base 

the  relationship  between  the  cradle  and  the  workpiece  for  a  constant  ratio 

of  roll 

U  generating  cone  surface  coordinate 

S  radial  location  of  cutting  head  in  coordinate  system 

r  radius  of  generating  cone  surface 

Equation  (2)  is  equivalent  to:  (3) 

Because  there  are  three  unknowns  U,  0,  and  4;  three  equations  must  be  developed  to  solve 
for  the  surface  coordinates  of  a  spiral  bevel  gear.  The  three  parameters  U,  0  and  <^<.  are 
defined  relative  to  the  cutting  head  and  cradle  coordinate  systems  (S^  and  S,)  respectively. 
These  parameters  can  be  transformed  through  a  series  of  coordinate  transformations  to  a 
coordinate  system  attached  to  the  workpiece.  Or  U,  0  and  can  be  mapped  into  X„,  Y„,  Z„ 


5 


in  coordinate  system  S*  attached  to  the  workpiece.  These  transformations,  used  in  conjunction 
with  two  other  geometric  requirements,  give  the  two  additional  equations. 

The  correct  U,  0  and  that  solves  the  equation  of  meshing,  must  also,  upon  transformation 
to  the  workpiece  coordinate  system  S,,,  result  in  a  axial  coordinate  Z„  that  matches  with  the 
value  of  Z  found  by  the  projection  of  the  tooth  in  the  XY  plane. 

-  Z  =  0  (4) 

This  equation  along  with  the  correct  coordinate  transformations  (see  Equation  11)  result  in  a 
second  equation  of  the  form: 

f^iU,  e,  <J)^  )  =0  (5) 

A  similar  requirement  for  the  radial  location  of  a  point  on  the  workpiece  results  in  the 
following: 

r  -  ^ {X^  +  Y^)  =  0  (®) 

This  is  shown  in  figure  2.1.  The  appropriate  coordinate  transformations  (see  Equation  11) 
will  convert  equation  6  into  a  function  of  U,  0  and 

f,iU,  e,  4)^)  =  0  (7) 

Equations  (3),  (5),  (7)  form  a  system  of  nonlinear  equations  necessary  to  define  a  point  on 
the  tooth  surface. 
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2.3  Solution  Technique 

The  three  equations  discussed  earlier  to  describe  the  tooth  surface  coordinates  are  nonlinear 
and  do  not  have  a  closed  form  solution.  They  are  solved  using  Newton’s  method  (ref.  10). 

An  initial  guess  Uo,  0o,  4>co  is  used  to  the  start  iterative  solution  procedure.  Newton’s 
method  is  used  to  determine  subsequent  values  of  the  updated  vector  (U^,  0^, 


yU 

= 

0J.-1 

+ 

yU 

yU 

(8) 


where  the  vector  Y  is  the  solution  of; 


3f2(er 

dU 

30 

3^(6^ 

du 

30 

30c 

3^3  (0*-^ 

3^3  (0f  ^ 

du 

30 

300 

■  ' 

►  =  < 

yk-1 

The  3x3  matrix  in  the  preceding  equation  is  the  Jacobian  matrix  and  must  be  inverted  each 
iteration  to  solve  for  the  Y  vector.  The  equation  of  meshing,  function  fj,  is  numerically 
differentiated  directly  to  find  the  terms  for  the  Jacobian  matrix.  Function  and  fa  cannot  be 
directly  differentiated  with  respect  to  U,  0  and  0,.  After  each  iteration  U''  *,  0'='*,  0^’'  *  (in  the 
cutting  head  coordinate  system)  are  transformed  into  the  workpiece  coordinate  system)  and  are 
transformed  into  the  workpiece  coordinate  system,  with  the  series  of  coordinate 
transformations  as  given  in  Equation  11. 
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(10) 


r  cotijr  -  U  cosilf 

V 

=  [W,„l 

U  sini|r  sin0 

.1. 

U  siinjj  cos6 

where, 

[M^^]  =[Ar^^/(<t)^)  [/Hp^]  [Af„^/((|)^)  ]  [M-SC] 

Each  matrix  [M]  above  represents  a  transformation  from  one  coordinate  system  to  another 
(ref.  4). 

Functions  fj  and  are  evaluated  by  starting  with  an  initial  U^,  0,,  and  performing  the 

transformations  in  Equation  (11)  and  evaluating  Equations  (4)  and  (6).  The  numerical 

differentiation  of  and  is  performed  by  transforming  Uk+inc,  G^+inc,  0k+*nc  (where  inc 
is  a  small  increment  appropriate  for  numerical  differentiation)  into  X^+inc,  Y^+inc,  Z„+iiic. 
Equations  (4)  and  (6)  are  then  used  to  evaluate  the  numerical  differentiation.  Function  f,,  fj, 
and  their  partial  derivatives  are  required  for  the  Jacobian  matrix  and  are  updated  each 
iteration.  The  iteration  procedure  continues  until  the  Y  vector  is  less  then  a  predetermined 
tolerance.  This  completes  the  solution  technique  for  a  single  point  on  the  spiral  bevel  gear 
tooth  surface.  In  this  way  the  coordinates  of  the  surface  of  the  tooth  are  defined.  This 
solution  technique  is  repeated  four  times  for  each  of  the  four  surfaces;  gear  convex,  gear 
concave,  pinion  convex  and  pinion  concave. 

Since  all  four  surfaces  are  generated  independently,  additional  matrix  transformations  are 
required  to  obtain  the  correct  tooth  thickness.  The  concave  surfaces  are  fixed  on  each  tooth. 
The  convex  surfaces  are  rotated  as  required.  The  angle  of  rotation  is  obtained  by  matching  the 
tooth  top  land  thickness  with  the  desired  value. 
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2.4  Gear  and  Pinion  Orientations  Required  for  Meshing 

After  generating  the  pinion  and  gear  surface  as  described  above,  the  pinion  cone  and  gear 
cone  apex  will  meet  at  the  same  point  (see  figure  2.2).  This  point  is  the  origin  of  the  fixed 
coordinate  system  attached  to  the  workpiece  being  generated.  To  place  the  gear  and  pinion  in 
mesh  with  each  other,  rotations  described  in  the  following  example  are  required  (ref.  5). 

1.  The  gear  tooth  surface  points  are  rotated  by  360/N,+  180  CW  about  the  global  Z,^,  axis, 
for  this  example,  the  rotation  is  190  degrees. 

2.  The  pinion  is  rotated  by  90  deg  CCW  about  the  global  Y  axis. 

3.  Because  the  four  surfaces  are  defined  independently,  their  orientation  is  random  with 
respect  to  meshing.  The  physical  location  of  the  gear  and  pinion  after  rotations  described 
above  correspond  to  the  gear  and  pinion  overlapping.  To  correct  this  condition  the 
pinion  is  rotated  CW  about  its  axis  of  rotation  until  surface  contact  occurs.  For  the 
example  used  in  this  study,  the  rotation  was  3.56  deg. 

To  make  a  complete  pinion,  the  generated  pinion  tooth  was  copied  and  rotated  12  times  and 
the  generated  gear  tooth  was  copied  and  rotated  36  times  with  the  aid  of  a  solid  modelling 
program  as  shown  in  figure  2.3. 
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CHAPTER  HI 


CONTACT  ANALYSIS  BY  THE  FINITE  ELEMENT  METHOD 


The  advantages  of  the  Finite  Element  Analysis  (FEA)  for  accurate  deformation  and  stress 
analysis  of  complex  forms  is  well  known.  This  chapter  provides  a  brief  outline  of  the  FEA 
analysis  carried  out  for  spiral  bevel  gears.  It  also  describes  the  fundamental  concepts  of  non¬ 
linearity  with  emphasis  on  automated  contact  analysis.  Since  the  research  reported  herein  is 
presented  using  the  general  purpose  finite  element  code  MARC  (ref.  11),  details  of  its  special 
features  used  for  the  analysis  and  the  description  of  various  blocks  of  the  input  deck  are  also 
discussed. 

3.1  FEA  of  Spiral  Bevel  Gears 

Only  a  few  researchers  have  investigated  finite  element  analysis  of  spiral  bevel  gears.  FEA 
analysis  has  been  done  on  a  single  spiral  bevel  tooth  using  an  assumed  contact  stress  distribution 
(ref.  4).  More  recent  FEA  spiral  bevel  gear  analysis  research  has  attempted  to  solve  the  contact 
stress  distribution  in  a  multi-tooth  model  (4  gear  and  3  pinion  teeth)  (ref.  3).  The  tooth  pair 
contact  zones  in  reference  3  were  modeled  with  gap  elements.  The  study  here  uses  software 
with  automated  contact  options  in  order  to  avoid  certain  limitations  in  the  use  of  gap  elements, 
such  as: 

(i)  It  is  difficult  to  connect  the  gap  elements  with  proper  orientations  across  the 
normal  from  one  surface  to  the  other  surface  parallel  to  the  contact  point. 

(ii)  The  orientation  of  the  contact  plane  remains  unchanged  during  deflection. 

(iii)  It  is  difficult  to  accurately  select  the  properties  such  as  appropriate  open/closed 
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stiffness  values,  selection  of  the  stiffness  matrix  update  strategy  and  efficient 
problem  restarts. 

New  advances  in  the  state  of  the  art  for  FEA  provide  deformable  body  against  deformable 
body  penetration  algorithms  which  can  be  used  to  establish  the  nonlinear  boundary  conditions 
for  contact  problems.  MARC  (ref.  11)  is  one  such  FEA  package  software  which  uses  this 
algorithm  to  automatically  detect  contact. 

3.2  Concepts  of  Nonlinear  Analysis 

In  linear  FEA,  a  simple  linear  relationship  exists  between  force  and  deflection  (Hooke’s 
law).  For  a  metallic  spring  under  small  strain,  the  force  F  equals  the  product  of  the  stiffness 
K  and  the  deflection  U. 

Also,  the  deflection  can  be  obtained  by  dividing  the  force  with  the  spring  stiffness.  This 
relationship  is  valid  as  long  as  the  spring  remains  linear  elastic,  and  the  deflections  are  such  that 
they  do  not  cause  the  spring  to  yield  or  break.  If  the  spring  material  is  changed,  for  example, 
from  steel  to  rubber,  the  linear  force-  displacement  relation  is  no  longer  valid.  It  becomes  a 
nonlinear  problem  i.e.  Ft^KU. 

A  nonlinear  structure  is  the  one  for  which  the  force-deflection  relationship  cannot  be  directly 
expressed  in  terms  of  a  set  of  linear  equations. 

The  three  major  types  of  nonlinearities  are: 

(i)  Geometric  nonlinearity  (large  deformations,  large  strains,  snap  through  buckling) 

(ii)  Material  nonlinearity  (plasticity,  creep,  viscoelasticity) 

(iii)  Boundary  nonlinearity  i.e.,  a  changing  status  (opening/closing  of  gaps,  contact, 
follower  force) 

A  nonlinear  system  cannot  be  analyzed  directly  with  a  linear  equation  solver.  However, 
it  can  be  analyzed  by  using  a  series  of  linear  approximations.  Each  linear  approximation 
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requires  a  separate  pass  or  iteration,  through  the  program’s  linear  equation  solver.  Each  new 
iteration  is  as  expensive,  in  terms  of  CPU  time,  as  a  single  linear  analysis  solution. 

In  the  preprocessing  phase  of  a  nonlinear  analysis,  everything  is  quite  similar  to  linear  FEA 
data  input  except  the  user  is  required  to  specify  certain  nonlinear  analysis  controls  (i.e.,  large 
displacements,  "contact"  parameters,  convergence  controls,  etc.)  and  additional  material 
properties  required  for  a  nonlinear  analysis. 

In  the  solution  phase,  of  nonlinear  FEA,  the  solver  must  perform  the  analysis  in  steps  or 
increments.  Within  each  increment,  the  code  will  also  iterate  as  necessary  until  equilibrium  is 
achieved,  before  proceeding  to  the  next  increment. 

In  the  post  processing  phase,  the  user  looks  at  quantities  like  stress  contours,  etc.  A 
nonlinear  analysis  takes  tens,  hundreds  and  sometimes,  thousands  of  increments,  thus,  usually 
requiring  a  high  speed  computer  with  lots  of  memory  for  a  reasonable  turnaround.  The 
objective  in  a  successful  nonlinear  analysis  is  to  obtain  a  converged  solution  at  a  reasonable  cost. 

In  large  deformation  analysis,  incremental  load  AF  and  displacement  AU  is  related  by  the 
tangent  stiffness  Kt-  In  solving  this  type  of  problem,  the  load  is  increased  in  small  increments, 
the  incremental  displacement  AU  is  found,  and  the  next  value  of  the  tangent  stiffness  is  found 
and  so  on.  A  brief  description  of  the  incremental  solutions  will  now  be  discussed. 

FEA  is  an  approximate  technique,  and  there  exist  many  methods  to  solve  the  basic  equations. 
In  nonlinear  FEA,  two  popular  incremental  methods  used  to  solve  the  nonlinear  equilibrium 
equations  are:  Full  Newton-Raphson  or  the  Modified  Newton-Raphson. 

The  Full  Newton-Raphson  Method  (see  figure  3.1)  assembles  and  solves  the  stiffness  matrix 
at  every  iteration,  and  is  thus  expensive  for  large  3-D  problems.  It  has  quadratic  convergence 
properties,  which  means  in  subsequent  iterations  the  relative  error  decreases  quadratically.  It 
gives  good  results  for  most  nonlinear  problems. 

The  Newton-Raphson  method  iterates  using  this  equation: 
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[Kj]  {At/}  =  {F^^}  -  {FJ 


(11) 


where, 

[Kj]  =  the  tangent  stiffness  matrix 

{AU}  =  the  displacement  increment 

{F«pp}  applied  nodal  force 

=  the  restoring  N-R  force  (the  loads  generated  by  the  elemental  stresses) 

=  the  residual  force 

The  program  updates  the  tangent  stiffness  matrix  [K-j-]  and  the  residual  Fjpp-F„  at  each 
iteration,  and  then  resolves  the  equation  given  above.  Convergence  is  achieved  once  F^pp-F^  is 
less  then  a  convergence  criterion  that  is  set.  If  F,pp  is  not  equal  to  F„,  the  system  is  not  in 
equilibrium. 

As  shown  in  figure  3.1,  the  1st  iteration  yields  a  displacement  AUj,  using  the  initial  elastic 
stiffness  and  the  applied  load  F,pp.  The  nonlinear  response  yields  a  force  value  F„  for  this 
displacement.  The  2nd  iteration  yields  AUj,  using  the  updated  tangent  matrix  and  the  residual 
load.  Subsequent  iterations  quickly  drive  the  analysis  to  a  convergent  solution.  This  solution 
guarantees  convergence  if,  and  only  if,  the  starting  AU  is  "near"  the  exact  solution.  This  could 

be  obtained  by  taking  smaller  load  increments. 

In  solving  the  convergence  problem,  one  cannot  neglect  the  general  FEA  conflict  of  expense 
versus  accuracy.  One  must  balance  computational  expense  against  accuracy.  Using  a  finer 
mesh  and  multiple  load  increments  can  often  lead  to  a  more  accurate  and  more  robust  (less  likely 
to  diverge)  solution,  but  usually  at  increased  expense. 

The  Modified  Newton-Raphson  method  does  not  reassemble  the  stiffness  matrix  during  every 
iteration  as  shown  in  figure  3.2.  It  costs  less  per  iteration,  but  the  number  of  iterations  may 
increase  substantially  over  that  of  the  full  Newton-Raphson  method.  It  is  effective  for  mildly 
nonlinear  problems.  In  this  analysis  the  full  Newton-Raphson  method  was  used. 
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3.3  Automated  Contact  Analysis 

Many  common  structural  features  exhibit  nonlinear  behavior  that  is  status-dependent.  The 
stiffness  of  these  systems  shifts  abruptly  between  different  values,  depending  on  the  overall 
status  of  the  item.  Status  changes  might  be  directly  related  to  load,  or  they  might  be 
determined  by  some  external  cause. 

Situations  in  which  contact  occurs  are  common  to  many  different  nonlinear  applications. 
Contact  forms  a  distinctive  and  important  subset  to  the  category  of  changing-status  nonlinearities 
(ref.  15).  Contact,  by  its  very  nature,  is  a  nonlinear  problem.  During  contact,  both  the  forces 
transmitted  across  the  area  and  the  area  of  contact  change.  The  force-displacement  relationship 
thus  become  nonlinear.  Usually,  the  transmitted  load  is  in  the  normal  direction.  In  this  report 
the  method  of  reference  1 1  is  used  to  perform  the  nonlinear  analysis. 

Reference  1 1  is  a  general  purpose  computer  program  for  linear  and  nonlinear  stress  and  heat 
transfer  analysis.  This  program  is  capable  of  solving  problems  with  nonlinearities  that  occur  due 
to  material  properties,  large  deformations,  or  boundary  conditions.  In  general,  the  solution  of 
nonlinear  FEA  problems  requires  incremental  solution  schemes  and  sometimes  requires  iterations 
within  each  load/time  increment  to  satisfy  equilibrium  conditions  at  the  end  of  each  step.  The 
program  features  relevant  to  gear  meshing  are  discussed  in  this  section. 

The  FEA  program  used  has  a  fully  automatic  CONTACT  option  which  enables  the  analysis 
between  finite  element  bodies  without  the  use  of  any  special  gap  or  contact  elements.  The 
procedure  is  capable  of  handling  the  following  types  of  contact; 

i)  deformable  bodies  against  rigid  surfaces 

ii)  deformable  bodies  against  deformable  bodies 

iii)  a  deformable  body  against  itself 

The  CONTACT  option  was  originally  designed  for  analysis  of  manufacturing  processes  such  as 
forging  or  sheet  metal  forming,  but  its  capabilities  have  been  expanded  to  meet  other  analysis 
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requirements.  The  work  presented  in  this  document  utilized  the  program  of  reference  11 
running  on  a  supercomputer. 

Contact  between  the  bodies  is  handled  by  imposing  non-penetration  constraints  (reference 
11).  The  non-penetration  constraint,  as  shown  in  figure  3.3  and  figure  3.4,  is 

Ua  •  n  <  D 

A  non-penetration  constraint  (that  the  surfaces  cannot  inter-penetrate)  is  usually  handled  in  FEA 
codes  by  one  of  three  techniques:  Lagrange  multipliers  (used  in  several  FEA  codes  that  offer 
"gap"  elements);  penalty  functions  (one  application  being  the  use  of  stiff  or  rigid  connecting 
members  to  approximate  the  constraint);  or  solver  constraints.  In  some  FEA  codes  which  offer 
explicit  dynamics  capabilities,  a  fourth  technique  is  the  direct  application  of  contact  forces.  Use 
of  a  "gap  element"  means  node  to  node  contact.  The  CONTACT  features  uses  the  solver 
constraints  approach. 

Solver  constraints  are  used  to  impose  the  non-penetration  constraint,  and  a  very  efficient 
surface  contact  algorithm  which  allows  the  user  to  simulate  general  2  and  3D  multibody  contact. 
Both  "deformable-to-rigid"  and  deformable-to-deformable"  contact  situations  are  allowed.  The 
user  needs  to  only  identify  which  bodies  might  come  in  contact  during  the  analysis.  Self¬ 
contact,  which  is  common  in  rubber  problems,  is  permitted.  The  bodies  can  be  either  rigid  or 
deformable,  and  the  algorithm  tracks  variable  contact  conditions  automatically.  Thus  the  user 
no  longer  needs  to  worry  about  the  location  and  open/close  status  checks  of  "gap  elements". 
Automatic,  in  this  context,  means  that  user  interaction  is  not  required  in  treating  multibody 
contact  and  friction,  and  the  program  has  automated  the  imposition  of  non-penetration 
constraints.  Also,  coupled  thermo-mechanical  contact  problems  (e.g  rolling,  casting,  forging) 
and  dynamic  contact  problems  can  be  handled. 

Real  world  contact  problems  between  rigid  and/or  deformable  bodies  are  actually  3D  in 
nature.  To  solve  such  contact  problems,  one  needs  to  define  bodies  and  their  boundary 
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surfaces.  The  definition  of  bodies  is  the  key  concept  in  analyzing  3D  contact  automatically. 
For  rigid  bodies,  one  can  define  such  surfaces  as:  4-point  patch;  ruled  surface;  plane;  tabulated 
cylinder;  surface  of  revolution;  and  Bezier  surfaces. 

Deformable  bodies  are  defined  by  the  elements  of  which  they  are  made.  Once  all  the 
boundary  nodes  for  a  deformable  body  are  determined,  4  point  patches  are  automatically  created, 
which  are  constantly  updated  with  the  body  deformation.  Contact  is  determined  between  a  node 
and  all  body  profiles,  deformable  or  rigid.  A  body  may  fold  upon  itself,  but  the  contact  will 
still  be  automatically  detected,  thus  preventing  self  penetration. 

The  user  must  define  bodies  so  that  their  boundary  surfaces  can  be  established.  Deformable 
bodies  are  defined  by  a  list  of  finite  elements,  and  rigid  bodies  are  defined  by  geometrical 
entities.  Because  the  contact  boundary  conditions  are  a  function  of  the  applied  load,  the 
analysis  must  be  carried  out  incrementally.  Within  each  load  or  time  increment  of  an  analysis, 
additional  iterations  may  be  required  to  stabilize  the  contact  zone.  Contact  problems  involve 
two  important  aspects: 

i)  the  opening  and  closing  of  the  gap  between  bodies 

ii)  friction  between  the  contacting  surfaces. 

The  MARC  program  establishes  a  hierarchy  between  the  bodies  so  that  at  a  given  contact 
interface,  one  body  is  the  contactor  and  the  other  body  is  the  contacted.  The  set  of  nodes  on 
the  boundary  surface  of  a  contactor  are  candidate  nodes  for  contact.  The  boundary  surface  of 
a  contacted  body  is  defined  by  a  set  of  geometrical  entities.  A  user  specified  contact  tolerance 
is  used  to  determine  the  body  separation  distance  which  determines  whether  two  bodies  are 
considered  to  be  in  contact  with  each  other.  The  contactor’s  boundary  nodes  are  prevented 
from  pienetrating  the  surface  of  the  contacted  body  by  imjxising  solver  constraints.  For  contact 
between  a  deformable  body  and  a  rigid  body,  a  displacement  constraint  is  applied.  For  contact 
between  two  deformable  bodies  MARC  applies  multipoint  constraints  in  the  form  of  ties.  The 
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ties  link  the  motion  of  one  node  on  the  contactor  body  to  two  adjacent  nodes  on  the  contacted 
body.  During  each  iteration  as  nodes  enter  and  leave  contact  or  slide  between  adjacent  node 
segments  on  contacted  bodies,  a  bandwidth  optimization  is  performed  to  reduce  the  computer 
processing  time  required. 

During  contact,  it  is  possible  to  move  bodies  around  during  an  analysis;  however,  the  user 
must  make  sure  that  deformable  bodies  do  not  have  rigid  body  modes.  A  minimum  number 
of  boundary  conditions  or  spring  element  attachments  must  be  applied  to  prevent  rigid  body 
motion. 

A  static  analysis  of  two  (or)  more  bodies  that  are  not  initially  connected  poses  special 
problems  with  a  finite  element  analysis,  if  one  of  the  bodies  has  a  rigid  body  motion  component. 
If,  at  any  time,  the  two  bodies  are  disconnected  then  the  stiffness  matrix  would  become  singular 
and  unsolvable  (in  a  static  analyses).  This  is  because  finite  elements  require  at  least  some 
stiffness  connecting  all  the  elements  together  along  with  sufficient  displacement  constraints  to 
prevent  rigid  body  motion.  In  order  to  overcome  this  difficulty,  weak  springs  have  been  added 
to  connect  the  bodies.  The  spring  stiffness  is  very  low.  This  stiffness  is  there  only  to  supply 
some  stiffness  to  the  system.  The  stiffness  should  be  negligible  compared  to  the  material 
stiffness,  so  that  it  has  no  effect  on  the  solution. 
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CHAPTER  rV 


FINITE  ELEMENT  MODEL  AND  ANALYSIS 


This  chapter  describes  the  procedure  for  assembling  a  finite  element  gear  pair  model  for 
analysis  with  reference  1 1  software.  Two  different  models  are  analyzed.  These  models  are 
generated  in  the  geometric  modelling  program  of  reference  14  (FEA  pre  and  post-processor). 
The  first  model  is  a  two  tooth  segment  and  the  second  model  is  a  seven  tooth  segment  of  a 
gearset.  Different  programs  used  to  create  these  FEA  models  and  the  boundary  conditions 
imposed  on  them  are  described.  Various  sets  of  rotations  carried  out  in  the  gearsets  for  the 
analysis  to  determine  the  contact  path  and  the  load  sharing  across  the  tooth  are  then  discussed. 
A  detailed  report  of  the  problems  faced  during  the  course  of  the  analysis  for  convergence  to  take 
place  and  memory  and  CPU  requirements  of  the  system  used  are  described.  Different 
parameters  of  the  finite  element  analysis  input  commands,  which  have  been  studied  to  reduce 
the  CPU  and  memory  requirements  of  the  system,  are  also  reported  in  this  chapter. 

4.1  Model  Descriptions 

To  model  spiral  bevel  gears,  the  machine  settings  and  the  gear  and  pinion  design  data  as 
given  in  Tables  I  and  II  are  necessary.  The  equation  of  meshing  and  the  kinematics  of  gear 
cutting  are  incorporated  in  the  computer  programs  to  generate  the  spiral  bevel  gear  model  in  the 
geometric  modeling  program  (ref  14).  The  input  data  (for  the  geometric  modeling  program) 
for  the  seven  tooth  model  is  obtained  from  references  16,  17.  This  is  a  10  x  10  mesh  input  to 
generate  4  gear  teeth  and  3  pinion  teeth  in  mesh  to  simulate  contact  on  multiple  teeth  of  a  spiral 
bevel  gearset.  The  input  also  includes  the  hub  attached  to  the  gear.  The  input  data  (for  the 
geometric  modeling  program)  for  the  two  tooth  model  is  obtained  from  reference  6  and  7.  The 
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programs  generate  a  NxM  mesh  on  the  tooth  profile  of  a  spiral  bevel  gearset. 


TABLE  I:  PINION  AND  GEAR  DESIGN  DATA 


PINION 

GEAR 

Number  of  teeth  pinion 

12 

36 

Dedundum  angle,  degree 

1.5666 

3.8833 

Addendum  angle,  degree 

3.8833 

1.5666 

Pitch  angle,  degree 

18.4333 

71.566 

Shaft  angle,  degree 

90.0 

90.0 

Mean  spiral  angle,  degree 

35.0 

35.0 

Face  width,  mm  (in) 

25.4(1.0) 

25.4  (1.0) 

Mean  cone  distance,  mm  (in) 

81.05  (3.19) 

81.05  (3.19) 

Inside  radius  of  blank,  mm  (in) 

5.3  (0.6094) 

3.0  (.3449) 

Top  land  thickness,  mm  (in) 

2.032  (0.080) 

2.489  (0.098) 

Clearance,  mm  (in) 

0.762  (0.030) 

0.92964  (0.0366) 

TABLE  II:  GENERATION  MACHINE  SETTINGS 


CONCAVE 

CONVEX 

CONCAVE 

CONVEX 

Radius  of  cutter,  r,  in 

2.9656 

3.0713 

3.0325 

2.9675 

Blade  angle,  ^|l,  degree 

161.9543 

24.33741 

58.0 

22.0 

Vector  sum,  L„ 

0.038499 

-0.051814 

0.0 

0.0 

Machine  offset,  E„ 

0.154576 

-0.1742616 

0.0 

0.0 

Cradle  to  cutter  distance,  s,  in 

2.947802 

2.8010495 

2.285995 

2.285995 

Cradle  angle,  q,  degree 

63.94 

53.926 

59.234203 

59.234203 

Ratio  of  roll,  M^., 

0,30838513 

0.32204285 

0.950864 

0.950864 

Initial  cutter  length,  u,  in 

9.59703 

7.42534 

8.12602 

7.89156 

Initial  cutter  orientation,  6,  degree 

126.83544 

124.43689 

223.9899 

234.9545 

Initial  cutter  orientation,  d>.,  degree 

-0.85813 

-11.38663 

-0.35063 

-12.3384 
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A  24  X  12  mesh  was  used  to  create  the  two  tooth  model.  Eight  noded,  isoparametric  HEX 
elements  were  used.  The  seven  tooth  model  and  the  two  tooth  model  are  as  shown  in  figures 

4.1  and  4.2.  The  seven  tooth  model  consisted  of  8793  elements  and  11261  nodes  and  the  two 
tooth  model  consisted  of  3116  elements  and  4452  nodes. 

4.2  Loading  and  Boundary  Conditions 

The  boundary  conditions  for  the  seven  teeth  and  two  teeth  models  are  as  shown  in  figures 
4.1  and  4.2.  The  boundary  conditions  are  applied  such  that  the  gear  teeth  are  made  to  pivot 
about  a  fixed  point  at  node  number  7872  for  the  seven  tooth  model  and  at  a  node  number  4448 
for  the  two  tooth  model.  In  both  models,  the  axis  of  rotation  for  the  gear  is  the  Z  axis.  The 
nodes  where  the  forces  are  applied  are  in  the  gear  hub  of  the  models.  Fixed  displacement 
boundary  conditions  are  applied  at  8  nodal  points  in  the  Z  direction  only  for  the  gear  and  in  all 
directions  for  the  pinion  as  shown  in  the  figures.  Since  this  is  a  static  problem  involving  two 
bodies  (the  pinion  and  the  gear)  in  contact,  as  described  in  a  previous  chapter,  weak  springs  are 
added  to  prevent  the  rotational  rigid  body  modes  for  the  gear.  Eight  springs  are  added  away 
from  the  contact  region.  The  springs  connect  the  comer  nodes  of  the  pinion  with  the  comer 
nodes  of  the  gear  on  the  faces  where  contact  occurs.  The  stiffness  of  the  springs  are  100  Ibs/in. 

This  is  insignificant  when  compared  to  the  tooth  contact  stiffness  and  therefore  does  not  effect 
the  overall  solution. 

The  maximum  torque  for  the  gear  mesh  studied  was  9508  in.»  lbs.  on  the  gear.  This  torque 
is  applied  as  a  concentrated  force  with  a  moment  arm  on  the  gear  hub.  This  concentrated  force 
for  the  seven  tooth  model  was  4714  lbs.  with  a  moment  arm  of  2.017  inches.  For  the  two  tooth 
model,  the  force  is  3392  lbs.  applied  at  a  radius  of  2.798  inches.  The  force  is  applied  in 
different  increments  for  convergence  to  occur.  The  details  of  the  incremental  loading  are 
discussed  later  in  this  chapter. 
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4.3  Model  Generation  to  Predict  the  Contact  Path 

To  predict  the  contact  path  of  the  spiral  bevel  gears  as  they  roll  in  and  out  of  mesh,  several 
rotations  have  been  carried  out  on  the  preliminary  model.  The  rotation  of  the  gear  and  the 
pinion  should  be  in  the  same  directions,  i.e.,  both  positive.*  For  the  model  being  analyzed, 
the  gear  has  36  teeth  and  the  pinion  has  12  teeth.  Therefore,  for  each  degree  of  gear  rotation 
the  pinion  has  to  be  rotate  three  degrees  to  avoid  interference  during  meshing.  For  the  seven 
tooth  model,  seven  such  rotations  were  carried  out.  For  the  two  tooth  model,  three  such 
rotations  have  been  carried  out.  From  the  positioning  of  the  preliminary  model,  the  rotations 
were  carried  out  in  both  positive  and  negative  directions  to  determine  the  contact  path.  In  the 
seven  tooth  model  the  pinion  was  rotated  by  +6, 0,-6, -12, -18, -24, -30  degrees  and  in  the  two 
tooth  model  the  pinion  was  rotated  by  -I- 6, 0,-6  degrees.  The  corresponding  rotation  of  the  gear 
were  -2,0, -f-2, +4, -1-6, +8,-1- 10  degrees  for  the  seven  tooth  model  and  -2,0, +2  degrees  for  the 
two  tooth  model.  (Zero  degrees  corresponds  to  the  initial  position  for  the  gearset  based  on 
solving  the  equation  of  meshing.) 

While  generating  the  models,  care  must  be  taken  not  to  change  the  node  numbers  and  the 
element  numbers.  This  eliminates  a  great  deal  of  editing  in  the  FEA  input  file.  The  command 
used  in  the  geometric  modeling  program  to  do  this  with  the  gear  and  pinion  part  names  is  as 
follows: 

Name,  (partname),  rot,  rotation  axis  data,  (partname) 

Care  must  also  be  taken  to  rotate  the  gear  about  the  Z  axis  and  pinion  about  the  X  axis 
directions  to  avoid  interference  of  the  gears  in  mesh. 


‘NOTE:  Although  all  gearsets  must  rotate  in  the  opposite  sense,  (i.e.,  if  the  gear  rotates  CW,  the  pinion  must 
rotate  CCW)  this  gearset  does  in  fact  rotate  in  the  same  direction  (i.e.,  both  gear  and  pinion  have  positive  rotation) 
based  on  the  right  hand  rule  and  the  defined  coordinate  system.  This  occurs  because  the  gear  rotates  about  the  Z 
axis,  with  the  positive  Z  axis  pointing  from  the  toe  of  the  gear  to  the  heel.  However,  the  pinion  rotates  about  the 
X  axis  with  the  positive  X  axis  pointing  from  the  heel  of  the  pinion  to  the  toe.  The  result  is  both  gear  and  pinion 
have  positive  rotation. 
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4.4  Assembling  the  Spiral  Bevel  Gear  Pair  for  FEA  Analysis 

The  procedure  to  assemble  a  spiral  bevel  gear  pair  to  perform  FEA  analysis  is  as  follows: 
STEP  1.  Initially  the  gear  model  is  generated  in  the  geometric  modelling  program,  from  the 
input  file  obtained  after  executing  the  computer  codes  from  references  6  and  7  as  discussed 
earlier.  The  model  is  optimized  and  the  node  and  the  element  id’s  are  compacted.  The  load 
and  the  boundary  conditions  are  applied  to  the  model. 

STEP  2.  A  neutral  file  is  created  for  translation.  This  creates  a  preliminary  data  input  file. 
This  data  file  needs  to  be  edited  for  the  non-linear  analysis.  Certain  commands  and  controls  are 
added  which  are  not  available  in  the  geometric  modelling  program.  These  changes  are 
discussed  in  detail  later  in  the  chapter. 

STEP  3.  After  editing  the  preliminary  data  input  file,  it  is  ready  to  be  submitted  for  the 
analysis.  (A  sample  input  file  is  given  in  appendix  A.  Appendix  B  gives  a  description  of 
specific  MARC  commands.)  To  submit  the  job  to  a  supercompter,  a  job  control  language  (JCL) 
file  needs  to  be  prepared.  (The  JCL  file  is  given  in  appendix  C.)  This  file  contains  the 
workspace  required  and  the  CPU  time  required  as  well  as  other  related  details  which  are 
discussed  later.  Also  a  user  subroutine  file  should  be  present  in  the  same  directory  as  the  input 
file.  This  file  suppresses  the  printing  of  input  in  the  output  file  generated  after  a  run. 

After  the  job  has  run,  four  files  are  created  which  are  as  follows: 

1.  List  file  (job.lst)  -  the  output  file  which  gives  listing  of  the  results. 

2.  Post  file  (job.post)  -  the  post  processing  file  is  used  as  input  to  the  geometric 
modeling  translator.  (MARPAT) 

3.  Restart  File  (job. restart)  -  written  if  this  option  is  included  in  the  input  deck. 
Useful  for  continuing  a  run  from  the  last  complete  run  or  any  complete  run  for 
which  the  post  data  is  asked  to  be  written. 
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4. 


Log  File  (user.O#)  -  gives  details  of  the  job  run  on  the  supercompter,  giving 
CPU  time  used  and  other  data.  It  also  gives  any  SYSTEM  errors  encountered 
while  the  job  is  running. 

If  there  are  any  errors  the  job  ends  and  the  list  file  exists  with  an  EXIT  number  depending 
on  the  error. 

STEP  4.  The  geometric  modelling  program  translator  is  used  to  convert  the  FEA  results  file 
into  readable  format.  The  post  file  is  translated  into  three  files  for  each  increment  of  the 
analysis  which  is  stored  in  the  post  file.  These  are  as  follows: 

1.  Element  File  (i#s0.els)  -  gives  the  elemental  stress  values. 

2.  Nodal  File  (i#s0.nod)  -  gives  the  nodal  stress  values. 

3.  Displacement  File  (i^sO.dis)  -  gives  the  displacement  values. 

4.  Log  files  (marpat.msg  and  marpat.crd) 

STEP  5.  To  view  the  results  in  the  postprocessor,  the  original  neutral  file  is  read  to  create  the 
model.  Then  results  are  checked  and  the  files  are  read  as  required  for  elemental,  nodal  or 
displacement  results. 

The  flow  diagram  to  show  the  steps  undergone  to  finally  prepare  the  model  for  the  analysis 
is  shown  in  figure  4.3.  A  typical  MARC  input  deck  is  shown  in  figure  4.4. 

4.5  Changes  Made  in  Preliminary  MARC  File  for  Nonlinear  Analysis 

Since  all  of  the  features  available  in  the  MARC  program  are  not  available  in  the  geometric 
modeling  program,  the  MARC  input  file  created  by  the  geometric  modeling  translator  must  be 
edited  to  add  the  following  commands  and  controls  for  the  gear  analysis: 

1.  SIZING  -  This  value  is  set  to  meet  the  greatest  workspace  requirement  for  a  given 
problem.  The  estimate  of  the  greatest  workspace  requirement  is  usually  based  on 
experience.  If  the  SIZING  value  is  not  enough  for  the  analysis,  the  FEA  program 
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automatically  switches  to  an  out  of  core  solution  mode.  The  out  of  core  solution  mode 
is  usually  not  preferred  and  should  be  avoided  as  out  of  core  disk  space  is  required  and 
CPU  time  increases.  Initially,  when  the  preliminary  FEA  input  file  is  obtained  SIZING 
has  a  value  of  400,000  words.  In  this  analysis  SIZING  is  changed  to  26  Mega  words 
(Mw)  for  the  two  tooth  model  and  to  28.5  Mw  for  the  seven  tooth  model.  With  this 
values  the  analysis  solutions  were  obtained  in-core.  Care  should  be  taken  in  specifying 
this  value.  This  value  should  be  1.7  Mw  less  than  that  specified  in  the  JCL  for  the 
memory  requirement.  This  accounts  for  the  loading  and  running  of  the  FEA  program. 
This  value  specified  could  be  1.7  Mw  less  than  the  JCL  value,  but  should  never  be  more 
than  the  JCL  value. 

2.  SETNAME  -  This  gives  the  number  of  items  in  defined  sets  for  elements  and  nodes. 
It  usually  allots  50  items  per  set.  Since  the  geometric  modeling  translator  converts  each 
named  component  into  a  set,  a  lot  of  sets  are  defined  with  elements  and  nodes  in  separate 
sets.  This  results  in  excessive  name  sets  being  defined.  These  defined  sets  are  not 
required  unless  it  is  necessary  to  define  the  nodes  and  elements,  which  are  difficult  to 
identify  from  the  model  but  are  easily  defined  in  these  named  sets.  In  that  case,  only  the 
desired  name  sets  are  kept  and  the  others  are  deleted.  This  not  only  reduces  the  input 
file  size  but  also  reduces  the  memory  requirements  to  run  the  job.  Because  with  an 
increase  in  name  sets,  the  workspace  memory  requirement  is  increased.  It  is  always 
better  to  reduce  the  workspace  required  to  permit  the  job  to  run  within  core.  Since  the 
models  used  have  numerous  named  components,  the  preliminary  FEA  input  file  has  many 
name  sets  with  element  and  node  sets  defined  separately.  The  SETNAME  has  a  value 
of  901  initially  which  is  changed  to  200  with  50  items  per  name  set  i.e.,  only  4  name 
sets  are  then  kept.  This  change  resulted  in  the  reduction  of  memory  requirements  for 
the  seven  tooth  model  run. 
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3.  DEFINE  SETS  -  These  are  kept  to  the  minimum  and  only  used  if  required  as  is 
discussed  in  SETNAME  change.  The  only  four  name  sets  GEARE,  PINIONE,  PROJ2E 
and  RIME,  which  define  the  elements  for  the  respective  regions,  are  used  for  defining 
individual  material  properties. 

4.  PRINT  -  This  parameter  is  specified  with  option  5  to  provide  output  of  nodal  contact 
information  as  nodes  contact  and  separate  from  surfaces  during  the  analysis. 

5.  POST  -  Post  processed  data  is  controlled  by  this  option.  The  type  of  post  process  data 
and  the  rate  at  which  it  is  recorded  is  specified.  The  post  codes  given  are  changed  to 
17,  131  and  133  which  represent  VON  MISES,  MAXIMUM  PRINCIPAL  and 
MINIMUM  PRINCIPAL  stresses.  The  number  of  increments  at  which  the  data  is 
recorded  is  added  in  column  45.  Also  0  in  column  20  is  changed  to  1  to  write  formatted 
post  data. 

6.  POINT  LOAD  -  This  option,  which  is  given  in  the  beginning  of  model  definition  cards, 
is  deleted  as  loads  specified  in  zero  (null)  increment  are  ignored. 

7.  CONNECTIVITY  -  In  column  15  to  suppress  element  connectivity  being  printed  in  the 
output  listing  a  numerical  number  1  is  added. 

8.  COORDINATES  -  In  column  20  to  suppress  the  nodal  coordinates  being  printed  in  the 
output  listing  a  numerical  number  1  is  added. 

Additional  MARC  commands  are  described  in  appendix  B. 
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4.6  Running  Large  Finite  Element  Analysis 

Initially,  many  trial  runs  were  submitted  to  study  how  each  of  the  options  behaved.  The 
variables  in  the  cards  for  each  option  were  changed  a  number  of  times  until  an  appropriate  one 
was  found.  Considerable  effort  was  spent  debugging  and  optimizing  the  finite  element  code 
running  on  the  supercomputer.  Because  this  problem  utilized  significant  computer  resources, 
much  specialized  computer  knowledge  was  needed.  This  section  attempts  to  document  some 
of  this  experience. 

In  the  OPTIMIZE  option,  the  method  used  was  changed  from  4,  the  Wave  front  method 
followed  by  Grooms,  to  use  option  2,  the  Cuthill  Mckee  method.  This  new  method  was  less 
expensive  and  saved  considerable  memory  requirements.  Also,  the  ELSTO  option  has  been 
added  for  this  analysis.  This  reduces  the  in-core  memory  requirements  below  the  28  Mega 
word  (Mw)  limit  defined  on  the  SIZING  card.  This  option  was  added  after  a  lot  of  memory 
shortage  requirements  were  faced.  Before  this  problem  was  resolved  more  memory  was 
required  than  what  could  be  given  in  the  SIZING  card  running  on  the  supercomputer.  The 
problem  finally  ran  with  much  less  than  200,000  blocks. 

Much  experimentation  was  done  on  the  CONTROL  option  variables.  This  option  sets  limits 
on  the  number  of  increments  and  recycles  which  may  be  performed  during  a  nonlinear  analysis. 

As  discussed  earlier,  spring  elements  are  used  to  prevent  rigid  body  modes.  The  SPRINGS 
option  is  a  list  of  nodes  and  their  degrees  of  freedom  which  are  connected  by  spring  elements. 
Nodes  were  identified  on  the  two  bodies  where  SPRINGS  were  to  be  connected.  A  number  of 
trial  runs  were  performed  to  check  the  effect  of  the  stiffness  value  given  to  the  springs. 

The  RESTART  option  is  used  to  save  and  retrieve  analysis  data  so  that  the  analysis  can  be 
restarted  for  additional  increments.  The  variables  were  adjusted  to  determine  how  the  loading 
could  be  continued  for  the  next  restart.  While  performing  the  initial  runs,  the  restart  option  was 
used  frequently.  This  was  because  the  number  of  increments  in  which  the  load  was  applied  to 
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reach  the  final  torque  in  the  gears  was  not  appropriate  for  the  solution  to  converge.  However, 
the  option  was  not  used  in  the  final  analysis,  since  an  efficient  the  number  of  load  increments 
was  eventually  found  for  efficient  convergence. 

The  loading  was  specified  incrementally  with  the  AUTO  LOAD  option.  AUTO  LOAD  will 
cause  the  current  load  vector  to  be  repeatedly  applied  (additively)  for  the  number  of  increments 
requested.  The  number  of  increments  to  be  specified  in  the  AUTO  LOAD  and  also  the 
CONTROL  card  sets  posed  another  problem.  The  actual  number  of  increments  to  be  run  were 
specified  in  the  AUTO  LOAD  card  set.  The  number  of  increments  in  the  CONTROL  card 
were  set  to  the  upper  limit  for  the  total  number  of  increments  allowed  for  an  analysis.  This  can 
be  left  blank  and  let  the  program  take  the  default  maximum  allowed.  The  POINT  LOAD 
option  is  used  to  apply  the  load  per  increment.  The  FEA  program  uses  incremental  values  of 
the  load  which  are  summed  to  give  the  total  value  that  is  used  in  that  increment.  The  last  value 
of  the  incremental  load  that  was  input  will  be  used  until  the  new  incremental  value  is  read  in 
to  replace  it. 

For  contact  analysis  with  AUTO  LOAD,  a  TIME  STEP  history  definition  card  set  must  be 
included.  Although  the  time  step  for  each  increment  is  arbitrary  for  this  research  analysis,  it 
must  be  included.  When  rigid  bodies  are  included  in  a  contact  analysis,  displacements, 
velocities  or  accelerations  are  specified  and  therefore  a  time  step  is  essential.  This  problem  has 
deformable  bodies,  not  rigid  bodies,  but  the  FEA  program  code  still  requires  a  TIME  STEP  card 
set.  An  arbitrary  time  step  of  10  units  per  increment  was  specified. 

In  the  CONTACT  option,  the  number  of  entities  present  in  body  1  and  body  2  have  to  be 
given  appropriately  depending  on  the  mesh  density  of  the  contacting  surface.  Care  has  to  be 
taken  to  define  the  two  bodies  as  body  1  and  body  2.  Identifying  elements  present  on  the 
contacting  surface  of  the  bodies  becomes  complicated  for  complex  models  like  the  spiral  bevel 
gears.  Many  times  during  this  research,  the  model  was  optimized  in  the  solid  modeler  due  to 
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some  small  modifications  on  the  model  which  caused  element  numbering  to  change.  This 
necessitated  the  element  numbers  to  be  identified  all  over  again. 

4.7  Convergence  Problems 

For  a  nonlinear  analysis,  many  load  steps  or  increments  must  be  used  to  reach  the  final  or 
fully  loaded  state.  In  the  process,  one  can  encounter  some  convergence  problems. 

A  torque  of  9508  in.  lbs.  was  applied  to  the  gear.  This  torque  was  input  as  a  series  of  point 
loads.  Since  this  is  a  nonlinear  problem,  the  entire  torque  should  not  be  applied  in  one 
increment.  Several  analysis  (about  60)  were  run  to  determine  the  sensitivity  of  the  structure 
meeting  the  convergence  criteria  as  a  function  of  the  size  of  the  applied  load. 

The  FEA  program  solves  the  non-linear  problem  on  an  incremental  basis  and  iterates  within 
each  increment  until  the  convergence  criteria  are  met.  These  criteria  are  specified  in  the 
CONTROL  data  block.  Too  many  iterations  within  an  increment,  is  an  indication  that  too 
much  loading  is  being  attempted  for  that  load  step.  Initially  the  gears  are  separated  by  a  small 
gap.  At  this  point  the  total  structure  has  very  little  stiffness.  Applying  any  torque  causes  the 
gear  to  rotate  and  contact  the  adjoining  pinion.  Very  small  load  increments  should  be  applied 
otherwise  the  matrix  updated  for  this  soft  structure  will  not  converge.  After  the  gear  contacts 
the  stiffness  of  the  structure  increases  drastically.  The  amount  of  loading  can  then  be  increased 
until  the  required  amount  of  torque  is  achieved.  For  the  seven  tooth  model,  a  total  of  8 
increments  were  used  to  reach  the  desired  torque,  while  for  the  two  tooth  model  a  total  of  10 
increments  were  used. 

While  performing  the  trial  analysis  some  increments  kept  repeating  as  the  convergence 
criteria  was  not  met.  Also  it  was  observed  that  load  was  automatically  decreased  by  the 
program.  The  conclusion  was  drawn  that  the  load  applied  was  too  high  for  the  solution  to 
converge,  which  was  seemed  correct  as  the  program  itself  was  decreasing  the  load  values.  This 


28 


gave  an  indication  that  the  load  step  needed  to  be  reduced. 

To  determine  the  correct  load  to  be  applied  many  trial  analysis  were  performed.  The  load 
was  reduced  from  about  100  lbs.  to  5  lbs.  but  the  solution  was  still  not  converging.  Each  time 
different  nodes  kept  coming  into  contact  and  disturbed  the  convergence  requirement.  The  nodes 
in  contact  were  found  to  lie  on  the  border  of  the  defined  contact  region.  Since  nodes  other  than 
those  defined  as  contact  region  nodes  were  trying  to  contact,  convergence  was  affected. 

The  next  step  taken  was  to  define  more  surface  nodes  on  the  contact  region  (CONTACT 
CARD).  After  this  the  convergence  criteria  was  changed  from  displacement  to  force  residuals 
(CONTROL  CARD).  This  eventually  led  to  the  convergence  of  the  analysis. 

A  summary  table  of  the  computer  runs  and  the  incremental  loads  given  to  the  seven  tooth 
model  and  two  tooth  model  is  given  in  Table  III. 


TABLE  III:  SUMMARY  OF  LOADING  AND  CPU  TIME 


SEVEN  TOOTH  MODEL 

Increment  Number 

Total  CPU  time  (sec) 

Load  increment  (lbs.) 

Total  Load  (lbs.) 

1 

437.5 

12 

12 

2 

1403.37 

2 

14 

3 

1968.24 

785 

799 

4 

2541.77 

785 

1584 

5 

3115.01 

785 

2369 

6 

4287.06 

785 

3939 

4886.40 


785 


4274 


TWO  TOOTH  MODEL 


II  -  , 


Increment  Number 

Total  CPU  time  (sec) 

Load  Increment  (lbs.) 

Total  Load  (lbs.) 

1 

121.11 

2 

2 

2 

150.38 

2 

4 

3 

181.13 

2 

6 

4 

212.98 

2 

8 

5 

309.7 

2 

10 

6 

405.28 

2 

12 

7 

775.31 

845 

857 

8 

856.49 

845 

1702 

9 

939.53 

845 

2547 

10 

1097.90 

845 

3392 

4.8  Supercomputer  Requirements 

Supercomputers  have  various  queues  with  different  memory  and  CPU  time  requirements  as 
shown  in  Table  IV.  Each  of  the  queues  has  either  1  or  2  jobs  running.  A  small  job  with 
SIZING  of  around  400,000  words  usually  takes  less  than  300  secs  of  CPU  time  and  can  be  in 
the  smallest  queue  called  debug.  But  as  the  size  of  a  problem  increases  the  CPU  requirements 
increase  which  necessitates  a  job  to  be  submitted  in  a  higher  queue.  The  2  jobs  analyzed  in  this 
research  required  26  Mw  and  28.5  Mw  of  memory  space  and  about  20  -  80  minutes  of  CPU 
time.  The  turnaround  from  each  job  was  very  slow.  A  job  on  an  average  would  take  a  day 
and  a  half  in  the  queue  before  running.  After  a  lot  of  experimentation,  the  job  was  first 
submitted  in  the  300  second  slot.  Because  of  this,  errors  if  any  would  show  up  in  10  minutes 
rather  than  in  2  days. 
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TABLE  IV:  TYPICAL  INPUT  QUEUES  FOR  SUPERCOMPUTER 


Day  Queue  Limits 

Night  Queue  Limits 

Weekend  Queue  Limits 

Queue 

Name 

ESH 

■sssm 

Max  Jobs  Runnine 

Max  Jobs  Running 

Max  Jobs  Running 

Flood 

Limit 

Class  User 

Class  User 

Class  User 

debue 

300 

30.2 

3 

1 

2 

1 

1 

1 

2 

1.200 

4.2 

2 

1 

2 

1 

1 

1 

2 

3.600 

4.2 

2 

1 

2 

1 

1 

1 

2 

7.200 

4.2 

2 

1 

2 

1 

1 

1 

1 

a‘ 

3.600 

8.2 

2 

1 

2 

1 

1 
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Day  Queues:  0800-1700  Mon.  thru  Thurs.  &  0800-2200  Fri. 
Night  Queues:  1700-0800  Mon.  thru  Fri. 

Weekend  Queues:  2200  Fri.  thru  0800  Mon. 

Allow  0.7  Mw  for  System  Memory  Overhead 


As  mentioned  earlier,  the  space  requested  in  the  job  control  data  (JCL)  set  should  be  1.7  Mw 
higher  than  that  given  in  the  SIZING  card.  This  is  required  for  the  FEA  program  to  be  loaded. 
A  detailed  listing  of  the  JCL  is  given  in  appendix  C. 

If  the  job  is  required  to  be  restarted,  the  input  and  output  RESTART  files  should  be 
mentioned  in  the  correct  format.  This  option  was  utilized  in  this  research  since  the  turnaround 
was  very  slow.  Only  a  few  increments  were  requested  for  the  analysis  and  were  saved  in  the 
restart  file.  By  reading  the  restart  file,  the  job  could  be  restarted  from  the  previous  increment. 

Computing  the  amount  of  workspace  required  by  the  FEA  program  is  a  complex  function 
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of  many  variables.  The  most  efficient  method  is  to  select  a  large  enough  workspace  to  handle 
a  variety  of  runs,  without  sacrificing  efficiency  or  wasting  core  space. 

Both  in-core  and  out-of-core  data  storage  schemes  are  available  in  the  FEA  program. 
Elements  may  also  be  stored  out-of-core  if  the  ELSTO  option  is  used.  The  FEA  program 
chooses  the  solution  type  automatically. 

The  in-core  solution  technique  is  used  when  the  workspace  size  specified  in  the  SIZING  card 
is  larger  than  the  total  workspace  needed  in  the  in-core  matrix.  When  the  workspace  required 
is  too  large,  program  uses  the  out-of-core  solution  techniques  and  show  how  much  space  each 
nodal  row  requires,  the  number  of  nodal  rows  per  buffer  and  how  large  each  auxiliary  file  would 
be.  The  buffer  size  can  be  increased  only  by  changing  allocation  on  the  SIZING  card.  The 
amount  of  size  to  be  given  is  based  on  experience.  For  large  problems,  the  exact  workspace 
requirements  can  be  determined  without  actually  running  the  job  by  inserting  a  STOP  parameter 
card  after  the  workspace  is  allocated. 

The  frequency  for  writing  restarts  and  POST  data  should  be  low  to  avoid  disk  space 
problems  since  these  files  are  very  large.  Initially  these  files  were  written  after  every  10 
iterations.  Finally  they  were  written  for  every  2  iterations  when  the  number  of  increments  to 
apply  the  load  was  reduced  for  the  job  to  run. 
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CHAPTER  V 


RESULTS 

Various  jobs  were  successfully  run  to  predict  the  contact  path  and  the  load  sharing  for  the 
double  tooth  contact  region  as  the  gears  roll  through  mesh. 

5.1  Elliptical  Stress  Distribution 

The  stress  distribution  in  the  contact  region  was  found  to  be  elliptical.  Figure  5. 1  shows  a 
typical  elliptical  contour  obtained  using  elemental  stress  values  on  the  pinion  surface.  Figure 

5.2  shows  typical  pinion  contact  stresses  with  nodal  values.  Because  of  the  large  nodal  stress 
gradient,  the  nodal  stress  ellipse  is  slightly  distorted  when  compared  to  the  ellipse  contour  with 
elemental  stresses.  Note  that  only  the  pinion  contours  and  the  stress  values  are  discussed  in  this 
research,  since  the  gear  teeth  share  a  similar  load  distribution.  Also  note  that  only  the 
minimum  principal  stresses  <u:e  recorded  at  the  contact  region  for  the  study. 

5.2  Gap  Element  and  Automated  Contact  Analysis  Comparison 

Contact  stresses  on  spiral  bevel  gears  were  studied  by  researchers  with  gap  elements  in 
references  16  and  17  on  a  similar  seven  tooth  model.  Comparison  of  the  results  with  automated 
contact  analysis  will  now  be  presented.  Both  models  contained  the  same  mesh  density, 
boundary  conditions,  material  properties  and  loading.  The  nodal  stress  results  of  pinion  tooth 
#1  obtained  from  gap  elements  and  automated  contact  analysis  are  as  shown  in  figures  5.3  and 
5.4.  Note  that  both  the  contours  show  the  highest  concentrated  stress  value  at  the  same  node. 
A  comparison  of  the  results  of  these  two  runs  are  as  follows: 
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Gao  Element 

Automated  Contact 
Analvsis 

Max  Nodal  stress 

-296,410  (psi) 

-291503  (psi) 

CPU  time  (approx) 

30  min 

80  min 

Elemental  stress 

-84,761  (psi) 

-113,577  (psi) 

Gap  elements  closed  or 
(nodal  points  with  contact) 

4 

8 

No.  of  increments 

4 

8 

The  number  of  contact  nodes  at  the  contact  region  in  the  automated  contact  FEA  program  was 
higher  than  identified  by  the  gap  element  FEA  program.  With  the  pinion  considered  body  1 
(contactor  body)  and  the  gear  considered  body  2  (contacted  body),  eight  nodes  contacted  as 
shown  above.  With  body  1  and  body  2  switched,  16  nodes  were  found  to  have  contact. 
Presumably  this  sort  of  discrepancy  occurred  because  the  mesh  was  too  coarse. 

5.3  Seven  Tooth  Model  Results 

As  discussed  earlier,  the  seven  tooth  model  pinion  was  rotated  from  +6  degrees  to  -30 
degrees  about  the  X  axis  and  the  gear  was  rotated  from  -1-2  to  -12  degrees  about  the  Z  axis. 
As  these  gears  were  rolled  there  was  a  shift  in  the  contact  region.  It  was  observed  that  as  one 
pinion  tooth  goes  out  of  mesh  and  the  load  on  it  reduces,  the  other  pinion  tooth  shares  more  load 
and  as  it  starts  coming  into  contact.  The  elliptical  stress  contours  for  pinion  tooth  #\  and  pinion 
tooth  #2  while  they  are  rotated  in  mesh  with  the  respective  gear  teeth  (for  all  the  rotations  stated 
above)  are  shown  in  figure  5.5. 

Figures  5.6  to  5. 19  show  plots  of  nodal  stresses  for  pinion  tooth  Ul  and  ffl  for  all  positions. 
Shift  in  the  contact  nodes  while  the  gearset  rolls  through  mesh  are  shown  in  figures  5.20  and 
5.21.  These  nodes  are  obtained  from  the  output  file.  The  contact  forces  on  the  pinion  and  the 
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gear  nodes  which  are  listed  at  the  end  of  the  last  increment  are  studied  to  identify  the  nodes  in 
contact.  The  contact  node  density,  or  all  the  nodes  that  contact  during  meshing,  is  as  shown 

in  figure  5.22. 

5.4  Two  Tooth  Model  Results 

The  two  tooth  model  was  generated  with  a  finer  mesh  with  24  nodes  along  the  length  of  the 
tooth  and  12  along  the  height  of  the  tooth. 

As  discussed  earlier,  the  two  tooth  pinion  is  rotated  in  six  degree  increments  to  the  following 
positions:  +6,  0,  and  -6  degrees  to  observe  the  shift  in  the  contact  eUipse.  The  direction  of 
rotations  were  similar  to  that  of  the  seven  tooth  model.  The  element  and  nodal  stress  results 
for  the  model  rotated  by  -6  degrees  are  shown  in  figures  5.23  and  5.24.  In  this  model  due  to 
the  flexible  hub,  the  gear  started  sliding  over  the  pinion  by  a  large  amount  resulting  in  more 
nodes  in  contact.  For  this  reason  the  contact  ellipse  was  sliding  towards  the  edge  of  the  pinion 
and  the  ellipse  contour  became  distorted.  A  summary  of  the  seven  tooth  and  two  tooth  models 
with  reference  to  modeling  and  analysis  are  given  in  Table  V. 


TABLE  V:  SUMMARY  OF  FEA  ANALYSIS 


FEATURES 

7  -  TOOTH  MODEL 

2  -  TOOTH  MODEL 

No.  of  elements 

8793 

3116 

No.  of  nodes 

11261 

4452 

No.  of  degrees  of  freedom 

33748 

13321 

"In-core"  Memory  required  (Mw) 

28.5 

26 

CPU  time  (seconds) 

4886.40 

1092.90 

No.  of  increments 

8 

10 

Contact  tolerance  (in.) 

0.0002 

0.0002 

Torque  applied  lbs.  in. 

9508 

9508 
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5.5  Contact  Path 


The  contact  path  is  defined  as  the  path  of  the  point  of  maximum  force  on  the  gear  tooth 
while  it  rolls  in  and  out  of  mesh.  These  values  are  determined  for  the  pinion  using  the  output 
listing  which  gives  the  contact  force  in  the  final  increment  before  converging.  Referring  to 
Figure  5.5,  it  should  be  noted  that  while  the  contact  ellipse  in  pinion  tooth  #2  shifts  from  the 
toe  of  the  tooth  to  the  center  of  the  tooth,  the  contact  ellipse  in  pinion  tooth  U\  shifts  from  the 
center  of  the  tooth  to  the  heel  of  the  tooth.  The  contact  path  for  pinion  tooth  #1  and  pinion  tooth 
ff2  are  plotted  in  figure  5.25.  The  two  curves  of  pinion  tooth  and  pinion  tooth  tf2  do  not 
overlap  because  the  mesh  density  is  too  coarse  and  also  the  rotations  are  not  small  enough. 

The  force  at  a  particular  node  is  taken  to  be  the  square  root  of  the  sum  of  the  squares  of 
the  forces  in  the  X,Y  and  Z  directions.  The  maximum  force  is  obtained  for  all  the  rotations 
of  the  model  and  are  plotted  on  the  pinion.  The  two  tooth  model  showed  a  shift  in  the  contact 
path  due  to  flexibility  effects.  The  hub  region  for  this  model  was  very  flexible  and  had 
deformed  a  lot  before  the  gear  pair  came  into  contact. 

5.6  Load  Sharing  in  Spiral  Bevel  Gears 

The  load  sharing  was  analyzed  in  the  double  tooth  contact  region  of  the  seven  tooth 
model.  The  total  forces  on  pinion  tooth  #1  and  pinion  tooth  #2  at  each  of  the  rotations  from  +6 
degrees  to  -30  degrees  for  the  seven  tooth  model  are  calculated.  This  is  done  using  the  contact 
table  obtained  from  the  output  listing  of  any  run.  The  total  contact  force  on  each  pinion  is 
determined  by 

/(EF,)*-  .  .  (£F,f  =  F„_  (“•) 

Table  VI  shows  the  load  on  pinion  tooth  ffl  and  pinion  tooth  if2  calculated  from  various  runs 
while  the  gears  were  rotated.  The  load  on  the  tooth  as  a  function  of  rotation  is  also  plotted  as 
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shown  in  figure  5.26. 


TABLE  VI:  MAGNITUDE  OF  CONTACT  FORCES 
IN  PINIONl  AND  PINION2  ACROSS  THE  CONTACT  REGION 


Degree  of  Rotation 

PINION  TOOTH  #1 

PINION  TOOTH  n 

of  Pinion 

F. 

n 

F. 

Fu 

F. 

m 

Fz 

Fa 

+6 

2870 

3294.4 

473.6 

4397.4 

0 

0 

0 

0 

0 

2710.4 

2759.0 

447.0 

3893.0 

154.8 

164.9 

29.9 

226 

-6 

2323.7 

2032.0 

391.17 

3110.98 

513.0 

515.0 

83.5 

767.7 

-12 

1655.0 

1908.4 

295.3 

2542.9 

1008.9 

1093.4 

164.5 

1495.91 

-18 

1114.0 

1267.1 

240.5 

1704.07 

1582.2 

1655.1 

175.5 

2296.16 

-24 

646.3 

763.2 

135.7 

848.2 

1873.3 

2019.4 

343.68 

2775.2 

-30 

200.91 

252.6 

45.8 

324.8 

2156.0 

2326.0 

450.4 

3203.33 
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5.7  Fillet  Stresses 


The  fillet  stresses  are  plotted  for  the  seven  tooth  model  with  load  sharing.  A  typical 
contour  of  these  stresses  for  one  particular  rotation  is  shown  in  figure  5.28.  Table  VII  gives  the 
maximum  principal  nodal  stresses  at  various  nodes  and  at  various  roll  angles.  The  location  of 
these  nodes  are  given  in  figure  5.27.  These  values  are  plotted  as  a  function  of  roll  angle  as 
shown  in  figure  5.29. 


TABLE  VII:  FILLET  STRESSES  AS  A  FUNCTION  OF  ROTATION 
FOR  DIFFERENT  NODES 


PINION  1 

Nodes— > 

883 

871 

856 

1745 

1733 

1718 

Rotations 

+6 

75966.13 

138493.3 

32738.5 

39224 

18491 

15248 

0 

44804.1 

109165.1 

37978.5 

45357.8 

16021.8 

12918.5 

-6 

24030.4 

66600.7 

29362.4 

67333.4 

17878.6 

10812.1 

-12 

10272.5 

21690.6 

12490.4 

74458.8 

28887.1 

9505.6 

-18 

3721.5 

1312.2 

5926.3 

81997 

55723 

11607.1 

-24 

414.4 

15953.0 

13148.4 

90294.8 

67987.5 

14639.2 

-30 

1806.3 

29486.5 

16799 

37915.3 

81093.5 

27288.9 
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CHAPTER  VI 


SUMMARY  AND  CONCLUSIONS 


A  prcxjedure  for  predicting  the  contact  path,  load  sharing,  contact  stresses  and  fillet  stresses 
of  spiral  bevel  gears  is  presented.  The  method  incorporates  the  following  steps; 

1.  A  model  was  developed  using  the  kinematics  of  the  manufacturing  process,  the  machine 
settings,  and  the  design  data  for  the  gearset.  The  model  was  generated  in  a  geometric 
modeling  program. 

2.  Automated  contact  analysis  option  in  a  nonlinear  finite  element  analysis  program  was 
used  to  perform  the  analysis. 

3.  Two  different  models  were  analyzed.  Model  I  consisted  of  three  gear  teeth  and  four 
pinion  teeth  in  mesh.  The  nodal  mesh  density  on  the  tooth  surface  was  10  by  10. 
Model  II  consisted  of  one  pinion  and  one  gear  tooth  in  mesh  with  a  nodal  tooth  surface 
mesh  density  of  24  by  12. 

4.  These  models  were  analyzed  by  rotating  the  gears  in  angular  increments  to  predict  the 
contact  path.  Model  I  was  rotated  a  total  of  36  degrees  and  Model  II  was  rotated  a  total 
of  12  degrees. 

5.  The  load  sharing  was  determined  from  the  contact  forces  on  the  nodes  at  the  end  of  the 
converged  solution  for  each  rotational  position. 

The  initial  FEA  results  for  Model  I  at  the  contact  region  compared  favorably  with  the  results 
by  earlier  researchers  using  gap  elements.  It  was  observed  that  in  the  double  contact  zone, 
when  the  contact  ellipse  in  the  (i)th  pinion  shifted  from  the  center  of  the  gear  tooth  to  the  heel 
of  the  gear  tooth,  the  contact  ellipse  in  the  (i+  l)th  pinion  shifted  from  the  toe  of  the  gear  tooth 
to  the  center.  The  load  for  each  tooth  was  calculated  as  the  sum  of  the  contact  forces  on 
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various  nodes  in  the  contact  region.  The  number  of  nodes  in  contact  changes  as  the  gear  rolls 
through  mesh. 

There  was  a  large  stress  gradient  between  adjacent  nodes  in  the  contact  region.  This 
indicates  a  need  for  mesh  refinement. 
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Figure  2.1  Projection  of  gear  teeth  over  XZ  plane 
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Figure  2 . 2  Gear  and  Pinion  Orientations  to  Mesh  with  Each  Other 


44 


Figure  2.3  3-D  Model  of  the  Spiral  Bevel  Gears 
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Figure  3.1  The  Full  Newton-Raphson  (N-R)  Method 
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Figure  3.2  The  Modified  Newton-Raphson  Method 
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BODY  1 


Contact  Tolerance 


Contacted 
Body  '  ! 


When  contact  is  detected,  Multi-Point 
Displacement  Constraint  Ties  are 
introduced  to  prevent  body  penetration 


Figure  3.4  Contact  Algorithm 
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Pinion  Tooth 


Gear  Tooth 


O  Degrees  of  freedom 
constrained 


Axis  of  Rotation 
for  Gear 


Figure  4.2  Two  tooth  model  with  boundary  conditions 
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User  Subroutines  for 
Analysis  Control 


Figure  4 . 3  Flow  diagram  for  the  MARC  analysis 
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PSI 


A  -6664. 
B  -18543. 
C  -30422. 
D  -42302. 
E  -54181. 
F  -66060. 
G  -77939. 
H  -89819. 
I  -101698. 
J  -113577. 
K  -125457. 
L  -137336. 
M  -149215. 
N  -161094. 
O  -172974. 


Figure  5.1  A  typical  elliptical  contour  in  the  pinion 

tooth  with  elemental  stresses  in  seven  tooth 
model  (PSI) 
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PSI 


A  85251. 
B  22459. 
C  -40334. 
D  -103126. 
E  -165918. 
F  -228711. 
G  -291503. 
H  -354295. 


Figure  5.2  A  typical  elliptical  contour  in  the  pinion  tooth 
with  nodal  stresses  in  seven  tooth  model 
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PSI 


-17629. 
-57455. 
-97281 . 
-137107. 
-176933. 
-216759. 
-256584. 
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5.3  Nodal  Stross  Result  on  Pinion  Obtained  Froni  the  Gap 
Element  Solution  and  Seven  Tooth  Model 
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Figure  5 . 4  Nodal  stress  result  on  pinion  obtained  from 

automated  contact  analysis  in  seven  tooth  model 
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Pinion  Tooth  #1 


-12  Degree  rotation 


Pinion  Tooth 


-18  Degree  rotation 


-24  Degree  rotation 


Figure  5.5  (continued) 
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A  109199. 


Figure  5.6  Nodal  stresses  on  pinion  tooth  #1  after  it  is 
rotated  by  +6  degrees  in  the  seven  tooth  model 
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Figure  5.7  Nodal  stresses  on  pinion  tooth  #2  after  it  is 
rotated  by  +6  degrees  in  the  seven  tooth  model 
(No  Contact) 
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Figure  5.8  Nodal  stresses  on  pinion  tooth  #1  after  it  is 
rotated  by  0  degrees  in  the  seven  tooth  model 
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Figure  5.9  Nodal  Stress  on  Pinion  Tooth  #2  after  it  is  rotated 
0  Degrees  in  the  Seven  Tooth  Model 
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Figure  5.10  Nodal  stresses  on  pinion  tooth  #1  after  it  is 
rotated  by  -6  degrees  in  the  seven  tooth  model 
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Figure  5.11  Nodal  Stress  on  Pinion  Tooth  #2  after  it  is  rotated 
-6  Degrees  in  the  Seven  Tooth  Model 
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Figure  5.12  Nodal  stresses  on  pinion  tooth  #1  after  it  is 

rotated  by  -12  degrees  in  the  seven  tooth  model 
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Figure  5.13  Nodal  Stress  on  Pinion  Tooth  #2  after  it  is  rotated 
-12  Degrees  in  the  Seven  Tooth  Model 
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Figure  5.14  Nodal  stresses  on  pinion  tooth  #1  after  it  is 

rotated  by  -18  degrees  in  the  seven  tooth  model 
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Figure  5.15  Nodal  Stress  on  Pinion  Tooth  #2  after  it  is  rotated 
-18  Degrees  in  the  Seven  Tooth  Model 
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Figure  5.16  Nodal  stresses  on  pinion  tooth  #1  after  it  is 

rotated  by  -24  degrees  in  the  seven  tooth  model 
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Figure  5.17  Nodal  Stress  on  Pinion  Tooth  #2  after  it  is  rotated 
-24  Degrees  in  the  Seven  Tooth  Model 
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Figure  5.18  Nodal  stresses  on  pinion  tooth  #1  after  it  is 

rotated  by  -30  degrees  in  the  seven  tooth  model 
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Figure  5.19 


Modal  Stress  on  Pinion  Tooth  «2  after  it  is  rotated 
-30  Degrees  in  the  Seven  Tooth  Model 
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Contact  Node  Density  for  the  Seven  Tooth  Model 
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Figure  5.23  Elemental  stress  results  for  two  tooth  model 
at  -6  degree  rotation 
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Figure  5.24  Nodal  stress  results  for  two  tooth  model  at 
-6  degree  rotation 
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Figure  5.25  Contact  Path  of  Seven  Tooth  Model  in  Pinion  Tooth  #1 
and  Pinion  Tooth  #2.  Shown  is  the  Location  of  the 
Highest  Contact  Force  For  Each  Rotation  Indicated. 
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Double  Contact  zone  Analyzed 


5.26  Load  on  Tooth  as  a  Function  of  Rotation 


Pinion  Tooth  #1 


Pinion  Tooth  #2 


Figure  5.27  Location  of  Nodes  Identified  in  Table  VII 


Pinion  Tooth  #1 


Figure  5.28  Typical  fillet  stress  contours. 

Shown  for  -6  degree  rotation  of  pinion 
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APPENDIX  A 


MARC  INPUT  DATA 


TITLE  BULK  DATA  CARDS  PRODUCED  BY  PATMAR  RELEASE  3 . lA 

TITLE  20-Jul-93  15:14 

TITLE  OPTIMIZED  7TEETH  MODEL 
S  •**  PARAMETER  CARDS  *♦• 

28500000  . . .  this  allocates  workspace  storage  memory 

200  ...  defines  max.  no  of  items  in  DEFINE  sets 

7  ...  element  type  is  8  noded  hexahedral 

gives  extra  information  about  contact  nodes 
. . .  reduces  in-core  memory  by  storing  elements  out-of-core 
signifies  end  of  parameter  cards 
S  ***  CONTROL  CARDS  FOR  POST-PROCESSOR  TAPE  •»* 

$  -  SPECIFY  POST  CODES  AND  LABELS  FOR  EACH  ELEMENT  VARIABLE. 

$...! - 1 - ! - 2 - : - 3 - : - 4 - ! - 5 _ ! _ 6 _ ! _ 7 _ ! _ 8 

POST 

$  1  in  20th  column  gives  formatted  post  file  which  is  readable  on  different  m/c's 
S  4  shows  the  frequency  at  which  POST  data  is  written 


SIZING 

SETNAME 

ELEMENTS 

PRINT, 5 

ELSTO 

END 


3  11  4 

17  VON  MISES  STRESS 

131  MAXIMUM  PRINCIPAL  STRESS 

133  MINIMUM  PRINCIPAL  STRESS 

j  *.*  model  DEFINITION  CARDS  *•• 

S  —  SPECIFY  DISPLACEMENTS  AND  THE  ASSOCIATED 


S. . . ! . 

_ 1. 

1 

_ 2... 

.  ! _ 3 

1 

_ 4 _ 

FIXED 

DISP 

3 

0.0 

3 

•  •  • 

2952 

2961 

2992 

3016  8500  8581 

9680 

9761 

0.0 

0.0 

0.0 

1 

2 

3 

.  ♦  • 

7872 

♦  .  . 

0.0 

0.0 

0.0 

1 

2 

3 

259 

263 

364 

380  2317  2333 

2528 

2561 

DOFS  AND  NODES  FOR  EACH  SET. 

.  .5 _ ! _ 8 _ ! _ 7 _ !  , 


.8 


three  sets  of  displacement  data 
prescribed  displacement 
Z  direction  constrained 
—  list  of  nodes  on  gear 

X/Y/Z  directions  constrained 
pivot  node  on  hub 


list  of  pinion  nodes 


$ 

$  —  A  NAMED  SET  MUST  BE  DEFINED  BEFORE  ITS  NAME  IS  USED  IN  THE  INPUT  DECK. 

$...! - 1 - ! - 2 - ! - 3 - ! - 4 - ! - 5 _ ! _ 6 _ ! _ 7 _ ! _ 8 

SElement  set  GE^ARE  is  th  list  of  elements  of  named  component  GEAR  in  PATRAN 

DEFINE  ELEMENT  SET  GEARE 


1837  1838  1839  1840  1841 
7426  7427  7428  7429  7430 
7441  7442  7443 

1842 

7431 

1843  1844 
7432  7433 

1845 

7434 

1846 

7435 

1847 

7436 

1848 

7437 

1849 

7438 

1850 

7439 

1851 

7440 

c 

c 

SDefines  PINION  elements 
DEFINE  ELEMENT  SET 

P INTONE 

1  2  3  4  5 

1831  1832  1833  1834  1835 

6 

1836 

7  8 

9 

10 

11 

12 

13 

14 

15 

c 

DEFINE  ELEMENT  SET 

PROJ2E 

5023  5024  5025  5026  5027 
6679  6680  6681  6682  6683 

5028 

6684 

5029  5030 
6685  6686 

5031 

6687 

5032 

5033 

5050 

5051 

5052 

5053 

c 

DEFINE  ELEMENT  SET 

RIME 

2116  2117  2118  2119  2120  2121 
8779  8780  8781  8782  8783  8784 
$  -  SPECIFY  CONNECTIVITY  FOR 

2122  2123  2124 
8785  8786  8787 
EACH  ELEMENT. 

2125 

8788 

2126 

8789 

2127 

8790 

2128 

8791 

2129 

8792 

2130 

8793 

c 

S...! - 1 - ! - 2 - ! - 3 _ ! _ 4 _ ! _ 5 _ ! _ 6 _ ! _ 7 _ ! _ 8 

CONNECTIVITY 


Sshows  total  no.  of  elements  in  the  model,  1  suppresses  connectivity  printout 
8793  1 


84 


3 

7 

54 

53 

3 

4 

59 

58 

8 

9 

4 

7 

55 

54 

4 

5 

60 

59 

9 

10 

5 

7 

57 

56 

6 

7 

62 

61 

11 

12 

8791  71124611248112581125211247112491126011256 

8792  71125111257112541125011255112611125911253 

8793  71125211258112571125111256112601126111255 
5  ___  SPECIFY  COORDINATES  FOR  EACH  NODE. 

S...! _ 1 _ ! _ 2 _ ! _ 3 _ ! _ 4 _ ! _ 5 _ ! _ 6 _ i _ 7 - ! - 8 

COORDINATES 

5#  of  directions  per  node,  total  no.  of  nodes,  suppression  of  coordinate  listing 
611261  1 

1  -3.429364  -0.538929  1.299387. 

2  -3.429364  -0.553653  1.292881 

3  -3.429364  -0.568376  1.286376 

4  -3.429364  -0.583099  1.279870 

5  -3.429364  -0.597823  1.273365 

11261  -1.094289  -0.899709  1.786157 

Sdefines  material  properties 
ISOTROPIC 

3  _  three  sets  of  material  properties 

1, VON  MISES 

3.0000E+07,0.3  _  defines  Young's  modulus,  poisson  ratio 

Slist  of  elements  not  nodes... 

SLIST  OF  GEAR  AND  PINION  ELEMENTS... 

GEARE  AND  PINIONE  . . .elements  defined  by  ELEMENT  SET  names 

2,  VON  MISES 
3.0000E+10,0.3 
$RIM  ELEMENTS.  .. 

RIME 

3,  VON  MISES 
6.0000E-t-13,0.3 

SPRJECTION  IN  FRONT  OF  THE  HUB.. 

PROJ2E 

SPRINGS 

258.1.2671.1.100  ...  Ist  node,  dof,  2nd  node,  dof,  stiffnes 

258.2.2671.2.100 

258.3.2671.3.100 

364.1.2626.1.100 

364.2.2626.2.100 

364.3.2626.3.100 

2528.1.5512.1.100 

2528.2.5512.2.100 

2528.3.5512.3.100 

2317.1.5541.1.100 

2317.2.5541.2.100 

2317.3.5541.3.100 


CONTACT 

2,450,450  ...no.  of  contacting  bodies  and  entities  on  their  surface 


,0.0002, , 


_ contact  tolerance  defined 


1,0 


_ definition  of  1st  body 


o;, 

Slist  of  gear  elements  not  nodes  on  the  contacting  surface . 

$GE_4  BACK  FACE  ELEMENTS . . 

2902  TO  2904  BY  1  AND  3109  TO  3111  BY  1  AND  3316  TO  3318  BY  4  AND  C 

3523  TO  3525  BY  1  AND  3730  TO  3732  BY  1  AND  C 

2666,2670,  AND  2673  TO  2685  BY  4  AND  2686  TO  2688  BY  1  AND  C 

2881  TO  2901  BY  4  AND  3088  TO  3108  BY  4  AND  C 

3295  TO  3315  BY  4  AND  3502  TO  3522  BY  4  AND  C 

3709  TO  3729  BY  4  AND  3745  TO  3771  BY  1  AND  C 
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SGE_: 

5  BACK  FACE 

ELEMENTS. .  . 

4180 

TO 

4182 

BY 

1 

AND 

4450 

TO 

4452 

BY 

1 

AND 

C 

4702 

TO 

4704 

BY 

1 

AND 

4945 

TO 

4947 

BY 

1 

AND 

c 

5215 

TO 

5217 

BY 

1 

AND 

5485 

TO 

5487 

BY 

1 

AND 

c 

4159 

TO 

4179 

BY 

4 

AND 

4429 

TO 

4449 

BY 

4 

AND 

c 

4681 

TO 

4701 

BY 

4 

AND 

4924 

TO 

4944 

BY 

4 

AND 

c 

5194 

TO 

5214 

BY 

4 

AND 

5464 

TO 

5484 

BY 

4 

AND 

c 

5734 

TO 

5739 

BY 

1 

AND 

5806 

TO 

5811 

BY 

1 

AND 

c 

5740 

TO 

5742 

BY 

1 

AND 

5815 

TO 

5820 

BY 

1 

2,0 

.defintion  of 

2nd 

body 

0., 

of  pi.ni,on  olemonos  not  nodos  on  the  contacting  surface . 

5PI_1  FRONT  FACE  ELEMENTS.. 

4  TO  16  BY  4  AND  40  TO  52  BY  4  AND  76  TO  88  BY  4  AND  C 

112  TO  124  BY  4  AND  148  TO  160  BY  4  AND  184,195,199,203,0 

282,291,300,303,362,369,376,383,460,465,470,475,0 
20,27,32,35,36,56,63,68,71,72,92,99,104,107,108,  AND  0 
128, 135, 140, 143, 144, 164, 171, 176, 179, 180, AND  0 
207, 214, 219, 222, 223, 306, 309, 314, 317, 318, AND  0 
385,387,389,392,393,  AND  480  TO  484  BY  1  AND  0 
SPI_2  FRONT  FACE  ELEMENTS... 

508  TO  520  BY  4  AND  580  TO  592  BY  4  AND  652  TO  664  BY  4  AND  0 

726  TO  738  BY  4  AND  796,807,811,815,903,912,915,988,995,0 

443,450,455,458,459,524,531,536,539,540,  AND  0 
596,603,608,611,612,668,675,680,683,684,  AND  0 
742,749,754,757,758,819,826,831,834,835,  AND  0 
918,921,928,929,930,997,999,1001,1004,1005,  AND  C 
1092  TO  1096  BY  1 
CONTACT  TABLE 

1  ...one  set  of  bodies  defined 

1  . . .bodyl  is  contacting 

2  ...body2  is  contacted 

OPTIMIZE,  2, ,,,  1  . . .Cuthill  Mcheee  bandwidth  optimizer  is  used 

20, 

CONTROL  ...sets  controls  for  iterations 

$1  of  increments,  max.  f  of  recycles,  min.  I  of  recycles,  residual  checking  invoked,. 
,  5, 0, 1, 0, , 

.15 

RESTART  ...optional  restart  command 

1,2 

PRINT  ELEM 

1,2 

STRAIN, STRESS 
1  TO  8793 

PRINT  NODE 

1,2 

LOAD , REAC , TOTA .STRESS 
SCONSTRAINT  POINTS... 

7872,2952,2961,2992,3016,8500,8581,9680, 9761,  C 
259,263,364,380,2317,2333,2528,2561,  AND  C 
SPI_1  FRONT  FACE  NODES . . . 

5  TO  20  BY  5  AND  55  TO  70  BY  5  AND  105  TO  120  BY  5  AND  C 
155  TO  170  BY  5  AND  205  TO  220  BY  5  AND  255,268,273,278,0 
384, 395, 406, 410, 493, 502, 511, 520, 632, 639, 646, 653, 751, AND  C 
25  TO  225  BY  50  AND  34  TO  234  BY  50  AND  C 

41  TO  241  BY  50  AND  46  TO  246  BY  50  AND  C 

49  TO  249  BY  50  AND  50  TO  250  BY  50  AND  C 

283,292,299,304,307,308,414,418,425,430,433,434,761,756,0 
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523,526,529,534,537,538,  AND  660  TO  670  BY  2  AND  C 
770  TO  780  BY  2  AND  C 
SPI_2  FRONT  FACE  NODES... 

585  TO  605  BY  5  AND  614,621,626,629,630,  AND  C 

705  TO  725  BY  5  AND  C 

805  TO  825  BY  5  AND  C 

905  TO  925  BY  5  AND  C 

734  TO  934  BY  100  AND  C 

741  TO  941  BY  100  AND  C 

746  TO  946  BY  100  AND  C 

749  TO  949  BY  100  AND  C 

750  TO  950  BY  100  AND  C 
1009  TO  1029  BY  5  AND  C 
1038,1045,1050,1053,1054,1117,  AND  C 

1130  TO  1145  BY  5, AND  1154,1161,1166,1169,1170,  AND  C 
1246, 1257, AND  1268  TO  1280  BY  4  AND  1287,1292,1295,1296,  AND  C 
1355,1364,1373,1382,1385,1388,1391,1396,1399,1400,  AND  C 
1494, 1501, 1508, 1515, AND  1522  TO  1528  BY  2  AND  1531,1532,  AND  C 
1613, 1618, AND  1623  TO  1638  BY  5  AND  1639  TO  1642  BY  1  AND  C 
$GE  4  BACK  FACE  NODES .  .  . 


3678 

TO 

3681 

BY 

1 

AND 

3948 

TO 

3951 

BY 

1 

AND 

c 

4198 

TO 

4201 

BY 

1 

AND 

4448 

TO 

4451 

BY 

1 

AND 

c 

4698 

TO 

4701 

BY 

1 

AND 

4948 

TO 

4951 

BY 

1 

AND 

c 

3652 

TO 

3677 

BY 

5 

AND 

3922 

TO 

3947 

BY 

5 

AND 

c 

4172 

TO 

4197 

BY 

5 

AND 

4422 

TO 

4447 

BY 

5 

AND 

c 

4672 

TO 

4697 

BY 

5 

AND 

4922 

TO 

4947 

BY 

5 

AND 

c 

4  972 

TO 

soil 

BY 

1 

AND 

C 

$GE  ! 

S  BACK  FACE 

NODES. 

.  .  . 

5488 

TO 

5491 

BY 

1 

AND 

5828 

TO 

5831 

BY 

1 

AND 

c 

6148 

TO 

6151 

BY 

1 

AND 

6448 

TO 

6451 

BY 

1 

AND 

c 

6778 

TO 

6781 

BY 

1 

AND 

7108 

TO 

7111 

BY 

1 

AND 

c 

5462 

TO 

5487 

BY 

5 

AND 

5802 

TO 

5827 

BY 

5 

AND 

c 

6122 

TO 

6147 

BY 

5 

AND 

6422 

TO 

6447 

BY 

5 

AND 

c 

6752 

TO 

6777 

BY 

5 

AND 

7082 

TO 

7107 

BY 

5 

AND 

c 

7412 

TO 

7421 

BY 

1 

AND 

7512 

TO 

7541 

BY 

1 

ERROR  ESTIMATE 

1,1 

END  OPTION  . . .  .teminates  model  definition  cards 

SHistory  or  Load  Definition  Cards . 

POINT  LOAD 


0.0,12.0,0.0 

5141 

AUTO  LOAD 
1 

TIME  STEP 

10.0 

CONTINUE 
POINT  LOAD 


0.0, 2. 0,0.0 

5141 

AUTO  LOAD 
1 

TIME  STEP 

10.0 

CONTINUE 
POINT  LOAD 


0.0,785.0,0.0 
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5141 

AUTO  LOAD 

6 

TIME  STEP 

10.0 

CONTINUE 

•  QP  INPUT  DECK' 


88 


APPENDIX  B 


Description  of  MARC  Commands 
The  MARC  input  deck  consists  of  three  blocks.  They  are  as  follows: 

Parameter  Cards: 

TITLE:  It  gives  the  title  given  to  the  problem.  This  problem  has  the  title  as  "7TEETH  MODEL". 

SIZING:  This  specifies  the  size  of  the  workspace  buffer  in  number  of  words.  This  size  should  be  1.7  Mwords  less 
than  that  specified  in  the  JCL  file.  This  is  required  for  MARC  program  to  run.  Here,  it  is  quite  a  large  value. 

ELEMENT:  This  gives  the  element  type  used  in  the  model.  The  one  used  here  is  number  7. 

PRINT, 5:  It  is  a  special  printing  option.  "5"  means  additional  contact  analysis  infonnation  regarding  nodes 
touching  or  separating  fiom  surfaces  is  given. 

END:  It  signals  the  end  of  Parameter  Card  block. 


Modd  Definition  Cards: 

These  cards  contain  the  FE  model  data  for  the  analysis.  The  data  represents:- 

a)  The  FE  mesh  topology  -  element  connectivity,  nodal  coordinates  and  sheet  thickness. 

b)  Material  properties 

c)  Loading  and  boundary  conditions 

d)  Nonlinear  analysis  control 

e)  Output  controls 

f)  Contact  analysis  controls 
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FE  mesh  topology 

It  is  defined  by  the  following  cards:- 

CONNECnVITY:  This  defines  connectivuty  for  the  elements  in  the  model.  A  typical  card  is  as  illustrated: 
1  7  12345678 

First  number  denotes  element  number. 

"7"  is  the  element  type.  It  is  HEX  8-noded  element. 

The  last  8  numbers  define  the  connectivity  of  the  element  in  the  counter  clockwise  direction. 

COORDINATES;  This  gives  the  coordinates  for  the  nodes  in  the  model.  A  typical  card  is: 

1  0.0  0.0  0.0 

"1"  means  1st  node  number 

The  other  three  numbers  give  the  X,  Y  and  Z  coordinates  of  that  particular  node. 

The  2nd  card  is 

6  in  1 

"6"  means  max.  number  of  coordinate  directions  to  be  read  in  per  node. 

"Ill"  means  the  total  number  nodes  in  the  model 

"1"  in  the  20th  column  suppresses  printout  of  the  nodal  coordinates  in  the  output  file. 

ISOTROPIC;  This  lets  one  define  material  prop>erties,  a  yield  criteria.etc. 

2nd  card: 

no.  of  sets  of  isotropic  material  data  to  follow. 

3rd  card: 
l.VON  MISES 

"1"  is  material  identification  set  no. 

"VON  MISES"  is  the  selected  yield  criteria 
4th  card: 

3.00E+07,0.3  gives  Young’s  Modulus  and  poisson  ratio 
Sth  card: 

1  TO  §  gives  list  of  elements  associated  with  this  material 

SPRINGS;  They  are  used  to  input  any  simple  linear  springs. 
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2nd  card: 


103,1,51,1,50 

"103”  gives  the  node  to  which  1st  end  of  spring  will  be  attached 
”1"  gives  the  DOF  at  above  node  to  which  spring  will  be  attached 
"51”  gives  node  to  which  the  other  end  of  the  spring  will  be  attached 
"1"  gives  DOF  at  above  node  to  which  spring  will  be  attached 
"SO"  gives  the  stiffiiess  of  the  spring. 

Loading  and  Boundary  Conditions 

FDfED  DISP:  It  is  used  to  prescribe  displacement  boundary  conditions. 

2nd  Card:  No.  of  sets  of  boudary  condition  cards  to  be  dehned. 

3rd  Card: 

"0.0, 0.0, 0.0”  give  the  prescribed  displacements  for  1st,  2nd  and  3rd  DOFs 
i.e  X,  Y  and  Z  directions. 

4th  Card; 

1,2 

"1*  means  the  X  direction  DOF 
"2*  means  the  Y  direction  DOF 
5th  Card: 

76,78,82  to  200  by  5.. 

This  gives  the  list  of  nodes  to  which  above  boundary  condition  are  applied. 

These  above  3  Cards  are  repeated  as  required  to  define  displacement  at  various  nodes. 


POINT  LOAD: 

2nd  Card:  no.  of  sets  of  point  loads  to  be  entered.  It  can  be  left  blank. 
Here,  it  is  left  blank. 

3rd  Card: 

0.0,10.0,0.0 

"0.0"  gives  nodal  load  associated  with  1st  DOF  i.e  X  direction 
”10.0"  gives  nodal  load  associated  with  2nd  DOF  i.e  Y  direction 
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"0.0”  gives  nodal  load  associated  with  3rd  DOF  i.e  Z  direction 
4th  Card;  Gives  list  of  nodes  having  the  point  load  given  above. 

*5141''  is  the  node  at  which  load  is  applied. 

Non-linear  Analysis  Control 

CONTROL:  It  lets  one  input  parameters  which  control  the  convergence  and  accuracy  of  the  non-linear  analysis. 
2nd  Card: 

204.0.1.0. 

*20"  means  max.  no.  of  load  steps 

”S’  means  max.  no.  of  recycles  during  an  increment 

"0"  min.  no.  of  recycles  during  an  increment 

’1*  this  flags  the  convergence  testing  on  displacements 

’0.0’  flags  for  relative  error  testing 

’blank’  default  Full  Newton  Raphson  iterative  scheme  is  used 
3rd  Card: 

’.15*  gives  a  relative  error  of  15% 

Output  Controls 

POST:  It  creates  a  post  processing  tape  for  PATRAN. 

2nd  Card: 

7  11  1 

’7’  is  no.  of  element  variable  to  be  written 

’1'  is  in  16-20  column  for  formatted  POST  TAPE 

’1’  in  20-25  column  is  to  write  connectivity  and  coordinates  on  POST  TAPE 

’1*  in  45  colunm  is  to  write  post  data  every  increment 

3rd  Card;  Gives  various  post  codes 

11-16  give  components  of  generalised  stresses 

17  gives  Equivalent  Von  Mises  stress 

PRINT  ELEM:  Gives  elements  at  which  output  is  to  be  printed 
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2nd  Caid: 


la 

*1'  is  no.  of  sets 

’2’  is  inclement  between  piintout.  Default  is  every  increment  printed. 

3id  Card: 

"STRAIN  STRESS"  are  values  to  be  printed 
4th  Card: 

1  TO  44  is  the  list  of  elements  to  be  printed 
Sth  Card:  it  lists  integration  points. 

PRINT  NODE:  gives  information  on  nodal  printout 
2nd  Card: 

1,2 

*1*  is  no.  of  sets 

*2*  is  increment  between  printout.  Default  is  every  increment  printed. 

3rd  Card: 

"TOTA, LOAD, STRESS*  are  values  to  be  printed 
4th  Card: 

SO  TO  85  is  the  list  of  nodes  to  be  printed 

ERROR  ESTIMATE:  It  gives  error  associated  with  FE  discretization.  Large  values  indicate  stress  gradients  are 
not  accurately  represented. 

2nd  Card: 

1,1 

"1"  in  1-5  column  is  for  stress  measure  to  be  evaluated 
"1"  in  6-10  column  is  for  geometric  measure  to  be  evaluated 

SUMMARY:  Gives  the  sumnuiry  of  ouQ)ut 

Contact  analysis  controls 
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CONTACT:  Allows  one  to  perfonn  automated  contact  analysis  without  use  of  GAP  elements  for  rigid  to 
deformable  contact  as  well  as  deformable  to  deformable  contact.  Here,  it  is  deformable-deformable  contact. 
2nd  Card: 

2,60,60 

”2*  tells  two  surfaces  (bodies)  will  be  defined 
”60'  shows  there  are  max.  of  60  entities  to  be  created  for  any  body 
"60”  shows  max.  no.  of  nodes  that  lie  on  the  deformable  surface 
3rd  Card: 

,0.001,  gives  the  distance  below  which  node  is  considered  touching  a  surface 
04th  Card; 

1.0 

”1”  1st  body 

”0"  implies  deformable  body 
Sth  Card  and  6th  Card: 
blank 
7th  Card: 

17  TO  32  gives  list  of  elements  for  body  one. 

Sth  TO  11th  Card  are  repeated  as  above  from  6th  Card 

CONTACT  TABLE;  It  is  used  for  deactivating  or  activating  bodies  when  the  CONTACT  option  is  used. 
2nd  Card: 

”1*  no.  of  sets  of  bodies  to  be  input 
3rd  Card: 

”1*  gives  the  touching  body  number 
4th  Card: 

”2”  list  of  bodies  for  which  the  above  body  will  detect  contact 

OPTIMIZE:  It  allows  a  choice  of  bandwidth  optimizers  to  be  invoked  .  Is  helpful  in  reducing  the  bandwidth 
1st  Card; 

OPnMIZE,4,„l 
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"4"  shows  that  Wavefront  algorithm  based  on  connectivity  followed  by  Grooms  algorithm  is  used.  It  is  an 
effective  technique  for  3-D  (complex)  meshes. 

2nd  Card: 

1  TO  44  gives  list  of  elements 
3id  Card: 

*3*  gives  acceptable  half-bandwidth  at  which  it  should  exit  down  the  iteration  loop. 

Load  Incrementation  Cards: 

TIME  STEP:  Allows  user  to  prescribe  time  step  for  static  analysis 
2nd  Card: 

"10*  means  time  step  of  10  sec 

AUTO  LOAD:  It  describes  number  of  equal  load  steps  applied 
2nd  Card: 

"2*  shows  2  equal  load  increments 

CONTINUE:  This  terminates  Load  Incrementation  or  History  Definition  Cards. 
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APPENDIX  C 


JCL  AND  USER  SUBROUTINE 


Th«  JCL  file 


#  OSER-tobibel  PW-tobibel 

#  QSUB  -r  kt7tr0 

#  QSUB  -IT  10800 

t  QSUB  -IM  30.2Mvr 
t  QSUB  -eo 
I  QSUB  -q  systems 
set  -X 

ja 

news 

dir"/hogs/tobibel/7 teeth 

job“7tr0 

cd  Sdic 

rm  Sdit/S job.lst  Sdir/S job. post 
Is  -alp 

marc“/wtk/vvmaEc/msrc)c5/marc 


requester's  name  could  be  different  from  user 
CPU  time  requested 
memory  size  requested 

written  only  for  priority  systems  queue 


qivas  news  about  MARC 

defines  where  the  input  file  is  present 
shows  the  job  name 

changes  to  the  directory  where  input  file  is  present 
deletes  old  output  and  post  files  before  beginning  new  run 

defines  the  directory  where  MARC  is  accessed  from 


$ 

SAccesses  MARC,  defines  input  file,  output  file,  and  post  file 
$If  RESTART  option  is  used,  then  output  restart  file  is  written 
Sand  also  for  subsequent  restart  the  input  restart  file  is  mentioned 
S 

Sraarc  i-Sjob.marc  o-$ job.lst  patran-S job. post  orestatt”$ job. restl  usub-Sjob.f  news-no 
Is  -alp 
ja  -at 


SUBROUTINE  INLIST 


MARC  USER  SUBROUTINE  TO  DISABLE  ECHO  PRINTING  OF  THE  INPUT  DATA. 


C 

WRITE  (6,6000) 
C 


RETURN 

6000  FORMAT  (  2(/),'  ENTRY  TO  -INLIST’  USER  SUBRODTINE  • ,  6X,  90  (' -  • ) , 

G  2(/),SX, 'ECHO  OF  THE  INPUT  DATA  HAS  BEEN  DISABLED’, 

G  2(/)/'  EXIT  FROM  "INLIST*  USER  SUBROUTINE SX, 90 (’»•), /  ) 

END 
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APPENDIX  F 

IMPORTANT  OUTPUT  USTING 
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key  to  stressj  strain  and  displacement  output 
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internal  element  variables 
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equal  load  incs  specified  for  1  increments 
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MARC-CRAY  K5-1,  OB/13/93,  output  for  increment  1.  with  2  prop. ids 
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